Effects of Diffusion Coefficients and Permanent Charge on Reversal Potentials in Ionic Channels

In this work, the dependence of reversal potentials and zero-current fluxes on diffusion coefficients are examined for ionic flows through membrane channels. The study is conducted for the setup of a simple structure defined by the profile of permanent charges with two mobile ion species, one positively charged (cation) and one negatively charged (anion). Numerical observations are obtained from analytical results established using geometric singular perturbation analysis of classical Poisson–Nernst–Planck models. For 1:1 ionic mixtures with arbitrary diffusion constants, Mofidi and Liu (arXiv:1909.01192) conducted a rigorous mathematical analysis and derived an equation for reversal potentials. We summarize and extend these results with numerical observations for biological relevant situations. The numerical investigations on profiles of the electrochemical potentials, ion concentrations, and electrical potential across ion channels are also presented for the zero-current case. Moreover, the dependence of current and fluxes on voltages and permanent charges is investigated. In the opinion of the authors, many results in the paper are not intuitive, and it is difficult, if not impossible, to reveal all cases without investigations of this type.


Introduction
Ion channels are proteins found in cell membranes that create openings in the membrane to allow cells to communicate with each other and with the outside to transform signals and to conduct tasks together [1,2]. They have an aqueous pore that becomes accessible to ions after a change in the protein structure that makes ion channels open. Ion channels permit the selective passage of charged ions formed from dissolved salts, including sodium, potassium, calcium, and chloride ions that carry electrical current in and out of the cell.
To unravel how ion channels operate, one needs to understand the physical structure of ion channels, which is defined by the channel shape and the spatial distribution of permanent and polarization charge. The shape of a typical ion channel is often approximated as a cylindrical-like domain with a non-uniform cross-sectional area. Within a large class of ion channels, amino acid side chains are distributed mainly over a "short" and "narrow" portion of the channel, with acidic side chains contributing permanent negative charges and basic side chains contributing permanent positive charges, analogous to the doping of semiconductor devices, e.g., bipolar PNP and NPN transistors.
The spatial distribution of side chains in a specific channel defines the permanent charge of the channel. The spatial distribution of permanent charge forms (most of) the electrical structure of the channel protein. The spatial distribution of mass forms the structure studied so successfully by molecular and structural biologists. Ions that move through channels are often only an Angstrom or so away from the permanent charges residing on acid and base side chains. In addition, electrical forces are in general much stronger than entropic forces. Thus, in most cases, the electrical structure is more important in determining how ions go through a channel than the mass structure. Sometimes, the dielectric properties ("polarization") of the protein contribute a charge that is significant. Then, the spatial distribution of dielectric properties becomes an important part of the electrical structure.
The most basic function of ion channels is to regulate the permeability of membranes for a given species of ions and to select the types of ions and to facilitate and modulate the diffusion of ions across cell membranes. At present, these permeation and selectivity properties of ion channels are usually determined from the current-voltage (I-V) relations measured experimentally [2,3]. Individual fluxes carry more information than the current, but it is expensive and challenging to measure them [4,5]. Indeed, the measurement of unidirectional fluxes by isotopic tracers allowed the early definition of channels and transporters and is a central subject in the history of membrane transport, as described in textbooks-for example, [6][7][8][9]. The precise definition and use of unidirectional fluxes are dealt with at length in the paper [5]. The I-V relation defines the function of the channel structure, namely the ionic transport through the channel. That transport depends on driving forces expressed mathematically as boundary conditions. The multi-scale feature of the problem with multiple physical parameters allows the system to have great flexibility and to exhibit vibrant phenomena/behaviors-a great advantage of "natural devices" [10]. On the other hand, the same multi-scale feature with multiple physical parameters presents an extremely challenging task for anyone to extract meaningful information from experimental data, also given the fact that the internal dynamics cannot be measured with present techniques. The general inverse problem is challenging, although specific inverse problems have been successfully solved with surprisingly little difficulty using standard methods and software packages [11].
To understand the importance of the relation of current and permanent charges, that is, the I-Q relation, we point out that the role of permanent charges in ionic channels is similar to the role of doping profiles in semiconductor devices. Semiconductor devices are similar to ionic channels in the way that they both use atomic-scale structures to control macroscopic flows from one reservoir to another. Ions move a lot like quasi-particles move in semiconductors. In a crude sense, holes and electrons are the cations and anions of semiconductors. Semiconductor technology depends on the control of migration and diffusion of quasi-particles of charge in transistors and integrated circuits. Doping is the process of adding impurities into intrinsic semiconductors to modulate its electrical, optical, and structural properties [12,13]. In a crude sense, doping provides the charges that acid and base side chains provide in a protein channel.
Ion channels are almost always passive and do not require a source of chemical energy (e.g., ATP hydrolysis) in order to operate. Instead, they allow ions to flow passively driven by a combination of the transmembrane electrical potential and the ion concentration gradient across the membrane. For other fixed physical quantities, the total current I = I(V, Q) depends on the transmembrane potential V and the permanent charge Q. For fixed Q, a reversal potential V = V rev (Q) is a transmembrane potential that produces zero current I(V rev (Q), Q) = 0. Similarly, for fixed transmembrane potential V, a reversal permanent charge Q = Q rev (V ) is a permanent charge that produces zero current I(V, Q rev (V )) = 0.
The Goldman-Hodgkin-Katz (GHK) equation for reversal potentials involving multiple ion species [14,15] is used to determine the reversal potential across ion channels. The GHK equation is an extension of the Nernst equation-the latter is for one ion species. The classical derivations were based on the incorrect assumption that the electric potential Φ(X) is linear in X-the coordinate along the length of the channel. This assumption is particularly unfortunate because it is the change in the shape of the electrical potential Φ(X) that is responsible for so many of the fascinating behaviors of transistors or ionic systems [16][17][18][19][20][21]. There was no substitute for GHK equations until authors of [22,23] recently offered equations derived from self-consistent Poisson-Nernst-Planck (PNP) systems, to the best of our knowledge.
In this work, focusing on basic understanding of possible effects of unequal diffusion coefficients and, as a starting point, we will use the classical PNP model with a piecewise constant permanent charge and a cylinder-like channel with variable cross-sectional area. The classical PNP model treats ions as point charges. Among many limitations, gating and selectivity cannot be captured by the simple classical PNP model. However, the basic finding on reversal potentials and their dependence on permanent charges and on ratios of diffusion constants seems important and some are non-intuitive and deserving of further investigation. In the future, more structural detail and more correlations between ions should be taken into considerations in PNP models such as those including various potentials for ion-to-ion interaction accounting for ion sizes effects and voids [24][25][26][27][28][29][30][31][32].
There have been great achievements in analyzing the PNP models for ionic flows through ion channels [5,28,[33][34][35][36], etc. Although mathematical analysis plays a powerful and unique role to explain mechanisms of observed biological phenomena and to discover new phenomena, numerical simulations are needed to fit actual experimental data and study cases where analytical solutions do not exist. Furthermore, numerical observations may give clues for more theoretical investigations. Indeed, numerical and analytical studies are linked; any progress in one catalyzes work in the other.
This paper is a mathematical study on some aspects of ionic flows via the PNP models. It uses established mathematical methods and analytical results [23,33] that are derived without further assumption from their underlying physical models. The numerical results, throughout the paper, are gained from the algebraic systems (15), (22), (23) and (27), obtained from reduced matching systems of analytical results in [23,33]. The nonlinear algebraic systems are then solved by the MATLAB R (Version 9.5) function fsolve that uses the trust-region dogleg algorithm. The trust-region algorithm is a subspace trust-region method and is based on the interior-reflective Newton method described in [37]. Our numerical results indicate that current-voltage and current-permanent charge and even zero-current relations depend on a rich interplay of boundary conditions and the channel geometry arising from the mathematical properties analyzed in [23,33,34,38]. Although the work here is presented in the context of biological ion channels, it is clear that the results apply to the artificial channels that are now being studied for their engineering applications.
The highlights of our studies in this paper as well as in [23,33,34,38] applied to the setup of this paper include: (i) a mathematically derived system for the zero-current condition (see System (15)) that can be used to determine the reversal potential in terms of other parameters (see Display (22)); (ii) an examination on how the reversal potential depends on permanent charge: its sign and its monotonicity in permanent charge (see Section 2.2); and a comparison between this reversal potential and that from GHK in the special setting (see Section 2.3); (iii) a characterization of monotonic dependencies of the reversal potential on the ratio of diffusion coefficients in terms of different conditions on the boundary concentrations (see Section 2.2), as well as effects of un-equal diffusion coefficients on signs of zero-current flux and its dependence on permanent charge (see Section 2.1); (iv) numerical spatial profiles under the zero current condition of the concentrations and electric potential, and hence the profiles of the electrochemical potentials for several choices of permanent charges that reveal special features of permanent charge effects (see Section 2.4, particularly, Remark 3); (v) numerical and analytical studies of I-V and I-Q relations, and zero-voltage current and its rich dependence on permanent charge (see Section 3.3).
Furthermore, there are several qualitatively important but non-intuitive results discussed in this work. These qualitative results may be helpful in guiding experimentation and some might not be apparent in intuitive thinking about ion channel behavior. Here are some examples: a. The zero-current flux J has the same sign as that of l − r (see Section 2.1). b. The magnitude of the ratio between of the two diffusion coefficients affects the monotonicity of the zero-current flux J in Q (see Section 2.1). c. I-Q curves are not monotonic in general (see Section 3.2). d. Rich phenomena of interplay between boundary conditions and diffusion coefficients in terms of monotonicity of zero-voltage current on permanent charge (Section 3.3).
To this end, we would like to emphasize that applying the geometric analysis allows us to identify and formulate quantities and properties that are crucial to biology, while also providing quantitative and qualitative understanding and predictions.
This paper is organized as follows. The classical PNP model for ionic flows is recalled in Section 1.1 to prepare the stage for investigations in later sections. In Section 2, we study zero current problems to investigate the corresponding fluxes and reversal potentials V rev . In particular, we compare a special case of the reversal potential with the GHK equation. Some other numerical observations are also provided to study profiles of relevant physical quantities in Section 2.4. In Section 3, we first recall the analytical results in [33] when diffusion constants are also involved. Then, numerical observations are provided to examine behaviors of current, voltage, and permanent charge with respect to each other in some general cases. Some concluding remarks are provided in Section 4.

Poisson-Nernst-Planck Models for Ionic Flows
The PNP system of equations has been analyzed mathematically to some extent, but the equations have been simulated and computed to a much larger extent [39][40][41][42][43]. One can see from these simulations that macroscopic reservoirs must be included in the mathematical formulation to describe the actual behavior of channels [24,44]. For an ionic mixture of n ion species, the PNP type model is, for k = 1, 2, ..., n, where − → X ∈ Ω with Ω being a three-dimensional cylindrical-like domain representing the channel of lengthL (nm =L × 10 −9 m), Q( − → X ) is the permanent charge density of the channel (with unit 1M = 1Molar = 1mol/L = 10 3 mol/m 3 ), ε r ( − → X ) is the relative dielectric coefficient (with unit 1), ε 0 ≈ 8.854 × 10 −12 F m −1 is the vacuum permittivity, e 0 ≈ 1.602 × 10 −19 C (coulomb) is the elementary charge, k B ≈ 1.381 × 10 −23 JK −1 is the Boltzmann constant, T is the absolute temperature (T ≈ 273.16 K =kelvin, for water); Φ is the electric potential (with the unit V = Volt = JC −1 ), and, for the k-th ion species, C k is the concentration (with unit M), z k is the valence (the number of charges per particle with unit 1), and µ k is the electrochemical potential (with unit J = CV) depending on electrical potential Φ and concentrations C k . The flux density − → J k ( − → X ) (with unit mol m −2 s −1 ) is the number of particles across each cross-section in per unit time, D k ( − → X ) is the diffusion coefficient (with unit m 2 /s), and n is the number of distinct types of ion species (with unit 1).
Ion channels have narrow cross-sections relative to their lengths. Therefore, three-dimensional PNP type models can be reduced to quasi-one-dimensional models. The authors of [45] first offered a reduced form, and, for a particular case, the reduction is precisely verified by the mathematical analysis of [46]. The quasi-one-dimensional steady-state PNP type is, for k = 1, 2, ..., n, where X is the coordinate along the channel, A(X) is the area of cross-section of the channel over location X, and J k (with unit mol s −1 ) is the total flux through the cross-section. Equipped with System (2), we impose the following boundary conditions, for k = 1, 2, · · · , n, One often uses the electroneutrality conditions on the boundary concentrations because the solutions are made from electroneutral solid salts, The electrochemical potential µ k (X) for the k-th ion species consists of the ideal component µ id k (X) and the excess component µ ex The excess electrochemical potential µ ex k (X) accounts for the finite size effect of ions. It is needed whenever concentrations exceed, say 50 mM, as they almost always do in technological and biological situations and often reach concentrations 1M or more. The classical PNP model only deals with the ideal component µ id k (X), which reflects the collision between ions and water molecules and ignores the size of ions; that is, where C 0 is a characteristic concentration of the problems, and one may consider 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 value problem (BVP) (2) and (3), then the electric current I is For an analysis of the BVP (2) and (3), we work on a dimensionless form. Set In terms of the new variables, the BVP (2) and (3) becomes, for k = 1, 2, · · · , n, with the boundary conditions Remark 1. The actual dimensional forms of quantities have been used for all figures throughout the paper, that is, and we take C 0 = 10 M,L = 2.5 nm and D 0 = 2.032 × 10 −9 m 2 /s, and, for diffusion constants [31],

Setup of the Problem
We now designate the case we will study in this paper. We will investigate a simple setup, the classical PNP model (9) with ideal electrochemical potential (5), and the boundary conditions (10). More precisely, we assume (A0) The ionic mixture consists of two ion species with valences where Q 0 is a constant.
We assume that ε > 0 in System (14) is small. The assumption is reasonable since, ifL = 2.5 nm = 2.5 × 10 −9 m and C 0 = 10 M, then ε ≈ 10 −3 [47]. The assumption that ε is small enables one to treat System (14) of the dimensionless problem as a singularly perturbed problem that can be analyzed by the theory of geometrical singular perturbations (GSP). GSP uses the modern invariant manifold theory from nonlinear dynamical system theory to study the entire structure, i.e., the phase space portrait of the dynamical system, and is not to be confused with the classical singular perturbation theory that uses, for example, matched asymptotic expansions.
We rewrite the classical PNP system (9) into a standard form of singularly perturbed systems and turn the boundary value problem to a connecting problem. We refer the readers to the papers [33] and [36] (with insignificantly altered notations) for details. Denote the derivative with respect to x by overdot and introduce u = εφ. System (9) becomes, for k = 1, 2, System (14) will be treated as a dynamical system with the phase space R 7 and the independent variable x is viewed as time for the dynamical system. A GSP framework for analyzing BVP of the classical PNP systems was developed first in [33,35] for ionic mixtures with two types of ion species. The model of ion channel properties involves coupled nonlinear differential equations. Until accomplished, it was not apparent that any analytical results could be found, let alone the powerful ones provided by geometrical singular perturbation. This GSP framework was extended to an arbitrary number of types of ion species successfully only when two special mathematical structures of the PNP system were revealed [36]. One special structure is a complete set of integrals (or conserved quantities) for the ε = 0 limit fast (or inner) system that allows a detailed analysis of a singular layer component of the full problem. It should pointed out that most of the integrals are NOT conserved for the physical problem since, no matter how small it is, ε is NOT zero. The GSP allows one to make conclusion about the BVP for ε > 0 small from information of ε = 0 limit systems. The other special structure is that a state-dependent scaling of the independent variable turns the nonlinear limit slow (or outer) system to a linear system with constant coefficients. The coefficients do depend on unknown fluxes to be determined as a part of the whole problem, and this is the mathematical reason for the rich dynamics of the problem. As a consequence of the framework, the existence, multiplicity, and spatial profiles of the singular orbits-zeroth order in ε approximations of the BVP-are reduced to a system of nonlinear algebraic equations that involves all relevant quantities altogether. This system of nonlinear algebraic equations defines the physical framework of the problem precisely. The system shows explicitly what has been guessed implicitly "everything interacts with everything else" and, in the cases analyzed in this paper, the system shows quantitatively how those interactions occur. This geometric framework with its extensions to include some of the effects of ion size [28,29,32] has produced a number of results that are central to ion channel properties [5,23,30,34,38,48]; for example, it was shown in [34] that a positive permanent charge may enhance anion flux as well as cation flux; and, in order to optimize effects of the permanent charge, the channel should have a short and narrow neck within which the permanent charge is confined; and, it was shown in [38] that large permanent charge is responsible for the declining phenomenon-decreasing flux with increasing transmembrane electrochemical potential. We refer the readers to the aforementioned papers for more details on geometric singular perturbation framework for PNP as well as concrete applications to ion channel problems.
In this paper, we will apply some results and follow the notations in [23,33] for analytical results where the quantities are all in their dimensionless forms. In addition, for simplicity, we use the letters l, r and Q 0 where l 1 = l 2 = l, r 1 = r 2 = r, Q 2 = 2Q 0 .

Remark 2.
We remind the readers that the quantities V, l, r, c k , Q, φ,μ k , J k , D k , and I are dimensionless quantities corresponding to the dimensional quantities V, L, R, C k , Q, Φ, µ k , J k , D k , and I, respectively, obtained from Display (8). We switch from dimensional form to the dimensionless form and vice versa several times throughout the paper. Dimensionless variables are convenient for illustrating and analyzing mathematical and general physical relations. Dimensional variables are necessary for showing how evolution has exploited those general relations.

Zero Current Problems with General Diffusion Constants
In this section, we study how boundary concentrations, electric potential, permanent charges, and diffusion constants work together to produce current reversal. Throughout this section, in order to express the effects of diffusion constants on zero-current flux and reversal potential, we study and compare the results for different cases of diffusion constants where D 1 = D 2 and where D 1 = D 2 , to indicate and emphasize the differences.
Diffusion is the phenomenon through which the spatial distribution of solute particles varies as a result of their potential energy. It is a spontaneous process that acts to eliminate differences in concentration and eventually leads a given mixture to a state of uniform composition. Fick's first law [49] describes diffusion of uncharged particles by ∂ t c = D∂ 2 xx c, where c is the concentration, D is the diffusion constant, and t is time. Frequently, the determination of diffusion constants involves measuring sets of simultaneous values of t, c, and x. These measured values are then applied to a solution of Fick's law to get the diffusion constants. Many techniques are available for the determination of diffusion constants of ions (charge particles) in aqueous solutions [31,[50][51][52][53], etc. When diffusion constants are equal, classical electrochemistry tells that many electrical phenomena "disappear" altogether, e.g., the "liquid junction" is zero. If the diffusion constants of potassium and chloride are equal, classical electrochemistry says that KCl acts nearly as an uncharged species. Indeed, this is the basis for the saturated KCl salt bridge used in a broad range of electrochemical experiments for many years. Therefore, the equal diffusion constants case is quite degenerate. Experimental measurements are exclusively performed under isothermal conditions to avoid deviation of D values. Nevertheless, even diffusion constants of certain ionic species may differ from one method to another, even when all other parameters are held constant. Everything becomes much more complicated mathematically when the diffusion constants are not equal, however. This complexity is what makes many biological and technological devices interesting, useful, and valuable. Some kinds of selectivity depend on the non-equality of diffusion constants as well.
Applying GSP theory to the classical PNP system (2) for two ion species with diffusion constants D k , k = 1, 2, the authors of [23] obtained an algebraic matching system with eleven equations and eleven unknowns for zero current problems and singular orbits on [0, 1]. They further reduced the matching system for the case where two ion valences satisfy z 1 = −z 2 . It follows that the reduced matching system for zero current where and, A is the geometric mean of concentrations at x = a, that is, and Note that, if h(x) is uniform, then H(x) is the ratio of the length with the cross-section area of the potion of the channel over [0, x]. The original of this quantity H(x) has its root in Ohm law for resistance of a uniform resistor. It turns out that the quantities α and β together with the value Q 0 are key characteristics for the shape and the permanent charge of the channel structure (see Section 4 in [34] for more detailed and concrete results about the roles of α and β on the fluxes).
To this end, we recall three relevant results from [23] on which most of our analytical and numerical studies are based.

Theorem 1 (Theorem 3.4 in
For fixed Q 0 and θ, the reversal potential V rev = V rev (Q 0 , θ) can also be determined and enjoy properties stated in the next two theorems. Recall that we denote J 1 = J 2 by J.
In what follows, numerical simulations are conducted with the help of analysis on System (15). The combination of numerics and analysis gives a better understanding of the zero-current problems and compliments some analytical results obtained in [23]. For our numerical simulations, we choose a = 1/3, b = 2/3 in Display (13) and h(x) = 1 for simplicity and for definiteness.

Zero-Current
We aim to clarify the relationships of ion fluxes with permanent charge and diffusion constants when current is zero.
Recall that fluxes J 1 and J 2 are equal for this case and let J denote it. For any permanent charge Q = 2Q 0 , once a solution (A, V) of System (15) is obtained, it follows from matching equations (see Appendix in [23]) that J is given by

Sign of Zero-Current Flux J
It was observed in [22] that the Nernst-Planck equation in Display (9) (with dimensionless quantities) gives, for k = 1, 2, Therefore, the sign of flux J k depends only on the boundary conditions l, r and V. Note that Equation (21) holds for any condition, not just zero-current condition.
For zero-current problem, V = V rev depends on l, r, D 1 , D 2 , and Q as well, in general. Thus, the sign of zero-current flux J seems to depend on all quantities and to be difficult to figure out. It is not the case. A consequence of Display (20) together with Theorem 1 is that: The zero-current f lux J has the same sign as that o f l − r The latter follows directly from Theorem 1 that, for zero-current, l − A has the same sign as that of l − r. This is consistent with observations in Figure 1 where D 1 = 1.334 × 10 −9 m 2 /s is fixed, and D 2 varies from the same value to D 2 = 2.032 × 10 −9 m 2 /s, and to a random large value. Concerning the dependence of the zero-current flux J on Q 0 , we have the following: (i) If D 1 = D 2 , then the zero-current flux J is an even function in Q 0 , and it is monotonic for Q 0 > 0.
In this case, θ = 0 and, hence, it follows from Theorem 1 that A is an even function in Q 0 and is monotonic in Q 0 for Q 0 > 0, and thus is the zero-current flux J from Display (20). (ii) If D 1 = D 2 , then the zero-current flux J is not an even function in Q 0 and the monotonicity of the zero-current flux J in Q 0 seems to be more complicated.
In this case, it can be seen that G 2 in Display (16) is not an even function in Q 0 , and, hence, the zero-current flux J is not. We would like to point out that, it follows from [38], for fixed ρ = D 2 /D 1 , no matter how large, one always has J → 0 as Q 0 → ±∞ that is consistent with the observations in Figure 1. (iii) Another fascinating result is that the magnitude of ρ = D 2 /D 1 affects the monotonicity of the zero-current flux J in Q 0 .
In this case, if one fixes D 1 , and let D 2 increases from small values to D 2 → ∞, (i.e., ρ → ∞), then it follows from Display (20) that there is a meaningful change in the monotonicity of the zero-current flux J, for small values of Q 0 that is not intuitive.
Let us consider the case where L < R and Q 0 < 0 is small. Recall that A is the geometric mean of concentrations at x = a. It follows from System (15) and (16) that, as Q 0 increases, (a) A increases if ρ ≈ 1 (that is θ ≈ 0), and consequently the zero-current flux J decreases; (b) A decreases if ρ 1 (that is θ 1), and, hence, the zero-current flux J increases.
Thus, depending on the size of ρ, the zero-current flux J may increase or decrease in Q 0 < 0, which is also consistent with the observations in Figure 1. The analysis for the case with L > R is similar.
It seems likely that the engineering, like evolution, will use these mathematical properties to control the qualitative properties of channels, technological, and biological.

Reversal Potential V rev .
Experimentalists have long identified reversal potential as an essential characteristic of ion channels [54,55]. Reversal potential is the potential at which the current reverses direction, i.e., V = Φ(0) − Φ(L) that produces zero current I. Using dimensionless form of quantities (see Remark 2), it follows from System (15) and (16) (where there are two ion species with valences z 1 = −z 2 = 1) that for general permanent charge Q = 2Q 0 = 0 with arbitrary diffusion constants [23], the variable A (the geometric mean of concentrations at x = a) can be solved uniquely from G 2 = 0 in System (15), and the reversal potential is then where B, S a , S b , and θ are defined in Displays (18) and (19).

Range of Reversal potential V rev
For fixed l, r, and for any given Q 0 , it follows from Theorem 2 that there exists a unique reversal potential V rev such that V rev ≤ | ln l r |. As Q 0 → ±∞, then V rev gets close to the boundary values, i.e., V rev → ± ln l r .

Zero Reversal Potential
One particular case is when the reversal potential is zero.
To examine under what conditions one obtains V rev = V rev (Q 0 ) = 0, it follows Theorem 2 that, Considering the second case above, the observations in Figure 2 show that, as ρ = D 2 /D 1 increases, magnitude of the corresponding Q 0 becomes larger. In fact, as ρ → ∞, then Q 0 → −∞.

Reversal Potential
For Q 0 = 0, one has V rev (0) = θ ln l r from Theorem 2. Therefore, then V rev (0) has the same sign as that of θ(l − r).
Let us consider the case where D 1 < D 2 . In that case, V rev (0) has the same sign as that of l − r. This is reasonable, since, for V = 0, we have |J 1 | < |J 2 | (since all but J k /D k are independent of D k in Equation (21)), and to help |J 1 | more than |J 2 | to get J 1 = J 2 for zero current conditions, one needs to increase V when l > r and decrease V when l < r , and that is why V rev (0) > 0 for l > r and V rev (0) < 0 for l < r . This is consistent with observations in Figure 2 as well. The analysis for the other case with D 1 > D 2 is similar.

It follows from Theorem 3 that
For θQ > 0, ∂ Q V rev has the same sign as that o f l − r.
This analytical result does not allow conclusions about the case for θQ < 0, however. The observations in Figures 2 and 3 show that the result holds for any θ and Q. Thus, we have Conjecture : V rev is increasing in Q f or l > r and decreasing in Q f or l < r.
We remark that, in Figure 3, we take L = 20 mM, R = 50 mM, and D 1 = 1.334 × 10 −9 m 2 /s and D 2 = 2.032 × 10 −9 m 2 /s which are diffusion constants of, say, Na + and Cl − , respectively (see the solid line), and D 1 = 1.334 × 10 −9 m 2 /s and D 2 = 0.792 × 10 −9 m 2 /s, where D 2 is the diffusion constants of Ca 2+ (see the dashed line). The reversal potential V rev is increasing in ρ i f l > r and is decreasing in ρ i f l < r.
This feature reveals a fantastic aspect that is not intuitive immediately. Recall Equation (21). Given the boundary values and diffusion constants, the values one obtains for all terms in Equation (21) except J k are independent of D k [36]. The relation surely holds for the zero-current condition, i.e., J 1 = J 2 with V = V rev . Now, let us fix D 1 and increase D 2 (so ρ is increasing). Then, |J 2 | increases since all but J 2 /D 2 in Equation (21) are independent of D 2 . Consequently, to meet the zero-current condition, we need to increase |J 1 |. Intuitively increasing V rev seems to lead to an increase in |J 1 |. This intuition agrees with the result for l > r. However, for the case with l < r, the result is the exact opposite of the intuitive result. That is, for l < r, it says, as ρ increases, V rev decreases. This counterintuitive behavior could be explained by the fact that c 1 (x) depends on V rev , and reducing V rev could still increase |J 1 |. In fact, l < r will result in reducing V rev , but c 1 (x) changes in a way that consequently increases |J 1 |.
To illustrate the counterintuitive behavior, we provide a numerical result in Figure 4. We choose C 0 , L and D 1 for Na + as mentioned in Remark 1. Now, suppose that D 1 2 = 0.792 × 10 −9 m 2 /s, and consider the boundary concentrations L = 20 mM, R = 50 mM and Q = 1 M. In this case, V rev = −16.7657 mV and J = −1.7632 × 10 −17 mol s −1 . Now, if we increase D 1 2 to D 2 2 = 2.032 × 10 −9 m 2 /s, which is Cl − diffusion constant, then V rev = −19.5527 mV and J = −1.8788 × 10 −17 mol s −1 . These values make sense now, based on the above discussion. Note that we just pictured the middle part of the channel in Figure 4 since the sides are almost identical. One should notice that it is hard to realize, from Figure 4, how L < R will result in reducing V rev . The complicated behavior discussed above convinces us that detailed analytical studies, even for special cases, could be critical for the design and interpretation of numerical results.  Figure 4. Graphs of C 1 (X) when D 1 is fixed, but we increase D 2 .

A Comparison with Goldman-Hodgkin-Katz Equation for V rev .
In this section, we first recall the GHK equation [14,15], which relates the reversal potential with the boundary concentrations and the permeabilities of the membrane to the ions. If the membrane is permeable to only one ion, then that ion's Nernst potential is the reversal potential at which the electrical and chemical driving forces balance. The GHK equation is a generalization of the Nernst equation in which the membrane is permeable to more than just one ion. The derivation of GHK equation assumes that the electric field across the lipid membrane is constant (or, equivalently, the electric potential φ(x) is linear in x in the PNP model). Under the assumption, the I-V (current-voltage) relation is given by For the case where n = 2 and z 1 = −z 2 = 1, the GHK equation for the reversal potential is The assumption that the electric potential φ(x) is linear is not correct when applied to channels in proteins. This is because proteins have specialized structure and spatial distributions of permanent charge (acid and base side chains) and polarization (polar and nonpolar side chains). Experimental manipulations of the structure of channel proteins show that these properties control the biological function of the channel. The GHK equation does not contain variables to describe any of these properties and so cannot account for the biological functions they control. A linear φ(x) is widely believed to make sense without channel structure presumably, in particular, where Q 0 = 0. However, this is not correct either. It follows from Formula (22) for Q 0 = 0 that the zeroth order in ε approximation of the reversal potential in this case is Figure 5 compares V rev (0, ρ) in Formula (24) with V GHK rev from the GHK-equation in Display (23). It can be seen from the left panel that, when l and r are not far away from each other (for example L = C 0 l = 20 mM, R = C 0 r = 50 mM), then the two curves have almost the same behavior. However, when we reduce L from 20 mM to 1 mM, then the right panel shows a significant difference between the two graphs.

Profiles of Relevant Physical Quantities
It follows that, for any given Q, once a solution (A, V) of Equations (15) and (16) is determined, all the other unknowns can be settled, and, hence, the approximation of the solution of boundary value problem can be obtained. We consider the dimensional form of quantities, and fix (Q, L, R, D 1 , D 2 ) to numerically investigate the behavior of C k (X) and Φ(X) throughout the channel. It follows from the Nernst-Planck equation that µ k (x) has the same sign as that of µ k (1) − µ k (0) or the opposite sign as that of J k ; in particular, µ k (x) is always monotonically increasing or decreasing. For the zero-current situation, the reversal potential depends on ALL other parameters; and so it seems that it would be hard to make general conclusions about µ k (x), for example, about its monotonicity. This is not true. In fact, in Section 2.1, we have concluded that the sign of zero-current flux J is the same as that of L − R, and, hence, µ k (x) has the opposite sign as that of L − R. For the case considered in this part, L = 20 mM < R = 50 mM, one has J < 0, independent of the value of Q. Therefore, µ k (x) > 0 for k = 1, 2, and, hence, µ k (x)'s are increasing for zero-current situation when L < R, independent of Q, as shown in Figure 8 for Q = 0.1 mM and in Figure 10 for Q = 1 mM. On the other hand, as one changes the value of Q, the profiles of concentrations c k (x)'s and electrical potential φ(x) may vary from monotone to non-monotone, as shown between Figure 7 for Q = 0.1 mM and Figure 9 for Q = 1 M.

Current-Voltage and Current-Permanent Charge Behaviors
Ionic movements across membranes lead to the generation of electrical currents. The current carried by ions can be examined through current-voltage relation or I-V curve. In such a case, voltage refers to the voltage across a membrane potential, and current is the flow of ions through channels in the membrane. Another important piece of data are current-permanent charge (I-Q) relation. Dependence of current on membrane potentials and permanent charge is investigated in this section for arbitrary values of diffusion constants.
To derive the I-V and I-Q relations, we rely on [33] where the authors showed that the set of nonlinear algebraic equations is equivalent to one nonlinear equation for A, the geometric mean of concentrations at x = a defined in Equation (17). All other quantities or variables such as fluxes, profiles of electric potential φ(x) and concentrations c k (x) can then be obtained in terms of A. It is crucial to realize that this is a specific result, not available for general cases. One can only imagine that the resulting simplification produces controllable and robust behavior that proved useful as evolution designed and refined protein channels. The reduction allowed by this composite variable can be postulated to be a "design principle" of channel construction, in technological (engineering) language, or an evolutionary adaptation, in biological language. In particular, the current I can be explicitly expressed in terms of boundary conditions, permanent charge, diffusion constants, and transmembrane potential in the special case that allows the determination of A. In what follows, we derive flux and current equations-when diffusion constants are involved as well-in terms of boundary concentrations, membrane potential, and permanent charge. The I-V, I-Q, J-V, and J-Q relations are investigated afterward in Section 3.2.

Reduced Flux and Current Equations
In this section, for simplicity, in addition to the assumptions at the beginning of the setup section (Section 1.2), we will also assume that h(x) = 1, a = 1/3 and b = 2/3, in particular, α = 1/3 and β = 2/3 (see Display (19)). It was shown in [33] that the BVP (9) and (10) can be reduced to the algebraic equation where B = l − A + r, S a , S b and N are defined in Display (18), and Once A is solved from Equation (25), we can obtain the flux densities and current equations as follows: For any given (l, r, D 1 , D 2 , Q 0 , V), there exists a solution for the flux J and current I. The numerical results in the next section give us more information on "current-voltage" and "current-permanent charge" relations.

Dependence of Current on Diffusion Constants
Now, we reveal a feature of the theoretical results that is not intuitive. Suppose that (l, r, Q 0 ) is given (V is still free and is allowed to take any value!). It follows from Display (18) for the definition of N that there exists an A so that N = 0. It consequently follows from Equation (26) , then η = 0. Therefore, from Display (27), For special values o f parameters (l, r, V, Q 0 ), the sign o f I is determined by the sign o f D 1 − D 2 . Figure 11 is a numerical simulation from Equations (25) and (26) of the I-V curves for several values of Q with D 1 = 1.334 × 10 −9 m 2 /s and D 2 = 2.032 × 10 −9 m 2 /s. One may suspect, based on the numerical observations, that the value of current I, obtained from Display (27), is unique for any V, and is monotonically increasing in V. However, this may not be correct, in general. This is important since the opening and closing properties of channels might be thought to arise from non-unique solutions [16,17].

I-V Curves and I-Q Curves
The sign of J k is determined by the boundary conditions, independent of the permanent charge. Nevertheless, as expected and seen in Figure 12, the magnitudes of J k 's, and, consequently, the sign and the size of the current I do depend on Q = 2Q 0 C 0 in general. (Here, Q would be the nonzero value of the permanent charge in dimensional form.) Treating V as a parameter, the current I is a function of Q. The numerical observations in Figure 12 indicate that, (i) there exists some V * such that, for V > V * , I(Q) has a unique maximum; (ii) there exists some V * such that, for V < V * , I(Q) has a unique minimum.
In particular, I-Q curves are not monotonic in general.  In addition, we claim based on numerical observations (not proven though) that there existŝ V (D 1 , D 2 ) = min{V * , |V * |}, such that (i) for any given V where |V | >V, the corresponding current I is non-monotonic in Q, but (ii) for any V where |V | <V, the corresponding current I is monotonic in Q.
(a) for D 1 D 2 + f (l,r) l+r > 0, the zero-potential current I l (ν) is decreasing in ν (increasing in Q) for l < r, and is increasing in ν (decreasing in Q) for l > r; (b) for D 1 D 2 + f (l,r) l+r < 0, the zero-potential current I l (ν) is increasing in ν (decreasing in Q) for l < r, and is decreasing in ν (increasing in Q) for l > r. Figure 13 illustrates some of the above conclusions. In addition, it suggests that the monotonicity of I(0) holds for all values of permanent charge, not only for small or large values. We emphasize that the monotonicity of current I with respect to permanent charge Q is just true for zero membrane potential, i.e., V = 0. Indeed, one should recall from Section 3.2 that, when V = 0, then the current I is not monotonic in Q.

Conclusions
In this paper, we first recall the analytical results in [23] for arbitrary diffusion constants. To investigate the reversal potential problem for which the current is zero, we do numerical investigations based on the analytical results in [23], where many cases are studied analytically. We derive several remarkable properties of biological significance, from the analysis of these governing equations that hardly seem intuitive.
Biophysicists are also interested in the relation of current-voltage (I-V), and current-permanent charge (I-Q), as well as reversal potential problems. To do that, we first recall the analytical results in [33], for arbitrary diffusion constants, to drive the flux densities and current equations explicitly. One way to characterize channels is the current at zero electric potential, that is, when V = 0, which has practical advantages. Since it is usually easier to measure a large current than a vanishing one, we analyzed this case as well. Furthermore, we briefly study the special cases of small and large permanent charge for zero voltage case, based on the analytical results of [34,38], respectively. To bridge between small and large values of permanent charges, we numerically study I-V and I-Q relations for this case as well.