**Physics of Ionic Conduction in Narrow Biological and Artificial Channels**

Editors

**Peter V E McClintock Dmitry G. Luchinsky**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Peter V E McClintock Department of Physics Lancaster University Lancaster United Kingdom

Dmitry G. Luchinsky Department of Physics Lancaster University Lancaster United Kingdom

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Entropy* (ISSN 1099-4300) (available at: https://www.mdpi.com/journal/entropy/special issues/ physics Ionchannels).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Entropy* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1646-2 (Hbk) ISBN 978-3-0365-1645-5 (PDF)**

The cover illustrates a partially hydrated potassium ion translocating through a sub-nanoscale pore in a monolayer graphene membrane. We thank Subin Sahu and Michael Zwolak for permission to use their image.

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**



## **About the Editors**

**Peter V E McClintock** is Emeritus Professor of Physics at Lancaster University. He was educated at Queen's University, Belfast (B.Sc.1962) and the University of Oxford (D.Phil. 1966). Following a postdoctoral appointment at Duke University, he came to Lancaster in 1968 where he has served in many roles including as Head of Department for six years. His contributions to the physics of ion channels, with collaborators, include the proposal of ionic Coulomb blockade as a mechanism of selectivity and its verification in models and experiments. His research interests also encompass the superfluidity of liquid helium, isotopic purification of He-4, stochastic dynamics of nonlinear oscillators, large fluctuations, rogue waves, wave turbulence, quantum turbulence, neutron dynamics in fluid helium, and the nonlinear dynamics of ageing, anaesthesia, autism, hypertension and other cardiovascular diseases. He is Editor-in-Chief of Fluctuation and Noise Letters.

**Dmitry G. Luchinsky** is Senior Research Fellow at Lancaster University and Senior Research Scientist with KBR Inc. He was educated at Moscow State University (M.Sc. 1983) and the Russian National Research Institute for Metrological Service (Ph.D. 1990). Following several appointments in the Institute, he came to Lancaster in 1995 and was transferred to the permanent staff in 2006. He is currently a contractor for KBR Inc at the NASA Ames Research Centre. He has been investigating the physics of ion channels for almost twenty years with contributions that include the statistical and linear response theory of ion channels and the energetics of ionic transitions through solid state nanopores. His research interests also include additive manufacturing in space and machine learning methods for the discovery of partial differential equations.

## **Preface to "Physics of Ionic Conduction in Narrow Biological and Artificial Channels"**

This is a book about ion channels. It has been written mostly by physical scientists and mathematicians, even though the most widespread and important manifestation of ion channels is in biology, where they are essential to life in all its forms. How do non-biologists get involved in such investigations? Everyone will have their own particular story but, for ourselves, it was the heady combination of scientific curiosity, a wish to contribute to the fundamental understanding of natural phenomena that clearly have crucially important applications, and the realisation that some of our physics knowledge and expertise might be relevant.

The key properties underlying ion channel function are those of conductivity and selectivity –the ability to select between different kinds of ions, allowing the favoured species to pass easily while blocking others. It is now appreciated that an understanding of this selective conduction requires physics, and that the physics of biological ion channels has a great deal in common with that of artificial nanopores. We came to realise that, in each case, there are intriguing analogies with the physics of electrons in quantum dots. Although the ability to predict the function of a channel from its structure remains elusive, recent advances have brought us tantalisingly close to a fundamental theory of ionic permeation, based on the statistical physics of ions within the channel.

The book provides a collection of cutting-edge papers on ionic permeation through narrow water-filled channels, both biological and artificial, reprinted from a recent Special Issue of Entropy. The invited authors were selected as being leading scientists in the field with whose work we were already familiar. In some cases they are our past or present collaborators, or collaborators of collaborators, in what is quite a specialised scientific area. The book describes the statistical physics of the permeation process, mathematical aspects, modelling by molecular dynamics, and experiments. The time is ripe for bringing together these mutually complementary approaches. We hope and believe that they will facilitate major breakthroughs in understanding, enabling the design of nanopores to meet particular technological requirements as well as improvements in drug design.

August 2021

**Peter V E McClintock, Dmitry G. Luchinsky** *Editors*

## *Editorial* **Introduction to the Physics of Ionic Conduction in Narrow Biological and Artificial Channels**

**Dmitry G. Luchinsky 1,2 and Peter V. E. McClintock 1,\***


"There is plenty of room at the bottom"

#### Richard Feynman

The permeation of ions through narrow water-filled channels is essential to life and of rapidly growing importance in technology. Reaching an understanding of the mechanisms underlying the permeation process requires an interdisciplinary approach, where ideas drawn from physics are of particular importance and have brought encouraging progress in recent years. This introduction sets into context the several ground-breaking papers presented in the *Entropy* Special Issue on "The Physics of Ionic Conduction in Narrow Biological and Artificial Channels".

Understanding, predicting and optimising the ionic selective transport properties of nanopores remains a critical challenge, both to nanotechnology and to biophysics. The last few decades have witnessed substantial progress in the analysis of such transport based on the use of a variety of experimental, numerical, and theoretical methods. Indeed, it would require several books to do full justice to the current state of the art in the field.

In some cases, the crystal structures (e.g., those of potassium, sodium, and calcium voltage-gated channels) have been discovered. This has provided invaluable insight, but has also thrown into sharp relief the structure–function problem: how to predict the conduction/selectivity properties of a known structure; or, conversely, how to design a structure with the required properties. A reliable solution to the problem promises to open new horizons in terms of pharmaceutical applications and the improved fabrication of solid-state nanopores for the sensing of molecules, desalination, DNA sequencing, and the other developments that together are marking a new era in nanotechnology.

Novel numerical methods and computer hardware now enable microsecond-long simulations of systems with hundreds thousands of atoms and the exploration of polarisable and quantum mechanical force fields. They provide unprecedented capabilities for reaching an understanding of experimental data and for the development of novel devices and techniques. Theoretical advances not only underlie many developments in molecular dynamics, including enhanced sampling and advanced force fields, but are also opening up new research frontiers and shedding fresh light on a number of longstanding problems such as binding probabilities, knock-on mechanisms of conduction, gating, electric double-layers, and local dielectric permittivity, just to mention a few.

It is now appreciated that selective conduction in biological ion channels has a great deal in common with that in artificial nanopores. In each case, there are intriguing analogies with the physics of quantum dots leading to the development of the theory of ionic Coulomb blockade. We dedicate this Issue to the memory of our late colleague, Dr Igor Kh. Kaufman, who developed an elegant theory of ionic Coulomb blockade in biological ion channels and suggested a simple classification of voltage-gated channels based on the charge of the selectivity filter.

At the same time, it is known that specific features of ionic conduction (e.g., dehydration, ion-specific binding affinities, protonation, the multicomponent and competitive nature

**Citation:** Luchinsky, D.G.; McClintock, P.V.E. Introduction to the Physics of Ionic Conduction in Narrow Biological and Artificial Channels. *Entropy* **2021**, *23*, 644. https://doi.org/10.3390/e23060644

Received: 14 May 2021 Accepted: 19 May 2021 Published: 21 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Department of Physics, Lancaster University, Lancaster LA1 4YB, UK; d.luchinsky@lancaster.ac.uk

of ion dynamics, the complex and adaptive structure of the ionic pathway, long-range interaction, local variation of the effective dielectric constant, highly correlated motion of more than one ion within a narrow channel, electric double layers, and water layering at the channel entrances) add many layers of complexity to the fundamental physics analogies.

This Special Issue brings together original high-quality papers on ionic permeation through narrow water-filled channels, both biological and artificial, from some of the best researchers in the field. It includes papers on the statistical physics of the process, on molecular dynamics and Brownian dynamics simulations, and on relevant experiments. Although any selection of papers can only be a narrow slice of the field, our aim is to emphasise the complexity and mutual interdependence of recent multifaceted progress in understanding the physics of ion channels and nanopores. The time is ripe for bringing together these complementary approaches, and we anticipate that they will facilitate major breakthroughs, enabling the design of nanopores to meet particular technological requirements as well as improvements in drug design and perhaps in personalised medicine.

Importantly, the Poisson–Nernst–Planck (PNP) and kinetic models remain among the principal tools for predicting current through nanopores, both in biology and nanotechnology. An example of the classical application of the PNP model to the analysis of reversal potentials and zero-current fluxes, in a system with a fixed profile of permanent charges and two mobile ion species, is provided by the paper by Mofidi et al. [1]. Rigorous analytic and numerical results establish the dependence of the electric and chemical potential profiles on voltage and permanent charge.

At the same time, it is well known that classical Poisson–Boltzmann (PB) and PNP theories do not take account of short-range ion–ion, ion–wall, or ion–water interactions in ion channels. Efforts to eliminate or ameliorate the effects of this deficiency of the continuum models have a long history. This stream of research is represented by the interesting paper of J.-L. Liu and R.S. Eisenberg [2], featuring the development of a molecular mean-field theory—a fourth-order Poisson–Nernst–Planck–Bikerman theory for modelling ionic and water flows in biological ion channels. The theory treats ions and water molecules, in channels of any volume or shape, with interstitial voids, polarisation of water, and ion–ion and ion–water correlations. It can be applied to electrolyte solutions in the nanopores of batteries and fuel cells.

The modelling of ionic currents with reduced models is extensively analysed by Boda et al. [3]. They show that channels are especially amenable to reduced modelling because their functions and the relationships between input parameters (e.g., applied voltage, bath concentrations) and output parameters (e.g., current, rectification, selectivity) are well-defined, allowing one to focus on the physics of input–output relationships rather than on the atomic-scale physics inside the pore. Based on decades of research, the authors propose four general rules for constructing good reduced models of ion channels and nanopores, focusing on the physics of input–output relationships rather than on atomic structure. The proposed rules relate to the importance of (1) the axial concentration profiles, (2) the pore charges, (3) choosing the right explicit degrees of freedom, and (4) creating the proper response functions. Examples demonstrating the application of these rules are provided. Further improvements in predicting the capabilities of reduced models can be achieved by incorporating into the solution of the one-dimensional electro-diffusion model the potential of the mean force obtained from MD simulations. The performance of two such methods is examined by A. Pohorille and M. A. Wilson1 [4] using stochastic simulations. These methods require neither knowledge of the diffusivity nor simulations at multiple voltages, which greatly reduces the computational effort needed to probe the electrophysiology of ion channels. They can be used to determine the free energy profiles from either forward or backward one-sided properties of ions in the channel, such as ion fluxes, density profiles, committor probabilities, or from their two-sided combination. In this work, large sets of stochastic trajectories were generated, individually designed to mimic the molecular dynamics crossing statistics of models for channels of trichotoxin, p7 from hepatitis C and a bacterial homolog of the pentameric ligand-gated ion channel

(LGIC). The authors found that the free energy profiles and the current–voltage curves obtained from the generated trajectories reproduce with good accuracy results obtained in molecular dynamics simulations.

The charged particles of which matter is composed move when an external electric field is applied, and their changed distribution is traditionally described in terms of a polarisation field. For insulators, it is usually possible to define a relative permittivity (dielectric constant) to quantify the material's responsiveness to the electric field. In ion channels, for example, the protein walls and the water are usually treated as dielectric continua with relative permittivities of around 2 and 80, respectively. This approach can be very helpful and revealing, but it involves greater approximation than spatial averaging because, as R.S. Eisenberg points out [5], the material's response to the electric field may be both nonlinear and time-dependent. In order to accommodate such phenomena, while simultaneously challenging physicists to review their knowledge of electromagnetism in biological dielectrics, he proposes and discusses an apparently minor change in Maxwell's first equation. It produces a major consequence when joined with Maxwell's second equation in that conservation of total current (including the displacement current) then emerges as a general principle. In one-dimensional systems like ion channels or electronic circuit components, the consequences are profound: there, total currents are equal at all locations at any given time, so the space variable does not appear in the description of total current.

There follow two papers reporting MD simulations of ion currents in biological and artificial channels. First, S.M. Cosseddu et al. [6] present an extended MD-based analysis of ion motion within the KcsA channel. They reveal complicated patterns of potassium currents that are governed by the structural variability of the selectivity filter. They show that ion motion involves the complex dynamics of a strongly correlated network of residues and water molecules. Intriguing features of self-organisation and readjustment of the network are analysed statistically and discussed in detail.

Secondly, we note that ionic transport in nano- to subnano-scale pores is highly dependent on translocation barriers and potential wells. These features in the free-energy landscape are primarily the result of ion dehydration and electrostatic interactions. For pores in atomically thin membranes, the ionic dynamics both inside and outside the geometrical volume of the pore can be critical in determining its transport properties. S. Sahu and M. Zwolak [7] examine regimes of transport that are highly sensitive to pore size due to the interplay of dehydration and interaction with pore charge, where picometer changes in the size (e.g., due to a minute strain), can lead to a large change in conductance.

We have already remarked upon the crucial importance of water, the electric doublelayer, water-layering, polarisation, and the resultant changes of local dielectric permittivity at the entrances of nanopores. Another approach to this problem is illustrated in the paper by T.-L. Horng [8]. Starting from the classical Helmholtz free energy functional for an electrolyte, including the solvation energies for anions and cations, the author follows the Bikerman modification by adding an entropy term to the functional, and he then extends the Bikerman approach by introducing ion-size-specific corrections to the theory.

The approach based on density functional theory (DFT), which works well near charged walls and in bulk electrolytes, can be extended to the analysis of the orientational ordering of water dipoles in membrane nanotubes. M. Drab et al. [9] analyse water ordering in nanotubes by minimising the corresponding Helmholtz free energy functional, also including the orientational entropy contribution of water dipoles, and deriving the modified Langevin–Poisson–Boltzmann (MLPB) model of the electric double-layer. The MLPB equation is solved in cylindrical coordinates to determine the spatial dependences of the electric potential, relative permittivity, and average orientations of water dipoles within charged tubes of different radii. Results show that for tubes of large radius, the macroscopic (net) volume charge density of co-ions and counterions is zero on the geometrical axis. This is attributed to effective electrolyte charge screening in the vicinity of the charged inner surface of the tube. For tubes of small radius, the screening region extends into the

whole inner space of the tube, leading to non-zero net volume charge density and non-zero orientational ordering of water dipoles near the axis.

The DFT results mentioned above are examples of statistical physics yielding insight into the function of ion channels and nanopores. This theme is continued and extended, first by Gibby et al. [10], who apply their recent derivation of an effective grand canonical ensemble and linear response theory of ion channels to analyse the conduction of the bacterial NaChBac selectivity filter. The authors compare their theory to experimental current–voltage and current–concentration dependences for a single channel and for a whole cell. They find that the statistical theory in the linear response regime correctly predicts many important properties of the NaChBac filter, including the concentration dependence of the reversal potential and the current–voltage relations. They also show that the theoretical results are consistent with MD simulations of the filter population at each binding site.

Secondly, the analysis of quantum mechanical effects in ion channels is another important direction, supported by the extended capabilities of modern quantum mechanics/molecular mechanics simulations. In this respect, interesting perspectives are opened by mapping the statistical mechanics of ion channels onto an effective quantum mechanics. Such investigations are reviewed by T. Gulden and A. Kamenev [11], who study the dynamics and thermodynamics of ion channels, considered as effective 1D Coulomb systems whose statistical mechanics is dominated by entropic effects that may be taken accurately into account by mapping onto an effective quantum mechanics. The corresponding semiclassical calculations for non-Hermitian Hamiltonians are conducted by applying tools from algebraic topology. The relationship of the solutions to the thermodynamics and correlation functions of multivalent solutions within long water-filled channels is discussed.

The actual properties of real nanopores are, of course, discovered by experiment, which has been leading the research in this area, especially since the discovery of the structure of the KcsA channel. In our Special Issue, experimental insight is provided by two of the leading research groups in the field.

O. Fedorenko et al. [12] discuss the properties of voltage-gated sodium channels (Navs). These channels play fundamental roles in eukaryotes but lack structural resolution, which renders understanding their structure–function relationships a challenging problem. Bacterial Navs, representing simplified homologues of their eukaryotic counterparts, have enabled both structural resolution and electrophysiological characterisation. However, their homotetrameric structure leads to an EEEE locus in the SF that is at odds with the DEKA locus of eukaryotic Navs. Indeed, prokaryotic Navs have long been considered more similar to eukaryotic calcium channels (Cavs) than to Navs, leading to the formulation of the "EEEE paradox". This was arguably solved by Kaufman et al. by the realisation that there is a critical D residue close to the EEEE ring of eukaryotic Cavs generating an effective EEEED locus of charge −5*e*. Fedorenko et al. present a follow-up of a previous study, aimed at mimicking the SF of eukaryotic Navs by engineering radial asymmetry into the SFs of bacterial channels. This goal was pursued with two approaches: co-expression of different monomers of the NaChBac bacterial channel in mammalian cells to induce the random assembly of heterotetramers, and the concatenation of four bacterial monomers to form a concatemer that can be targeted by mutagenesis on specific strands of the SF, thereby introducing asymmetry. Patch-clamp measurements and MD simulations showed that an additional gating charge in the SF leads to a significant increase of Na<sup>+</sup> and a modest increase in Ca2<sup>+</sup> conductance in the NavMs concatemer, in agreement with the behaviour of the population of random heterotetramers with the highest proportion of channels with charge −5*e*. This study confirms that, although the charge at the SF is important, it is not the only factor affecting conduction and selectivity. It also offers new tools extending the use of bacterial channels as models of eukaryotic ones.

The work by A. Chernev et al. [13] reviews the most promising approaches to the fabrication of artificial nanofluidic devices capable of reproducing the properties of single ion channels. It is shown that modern technologies have great potential in allowing one to test various theoretical models of ion channels. The review aims to highlight ionic Coulomb blockade—the phenomenon which (see above) can often be a key player in ion channel selectivity. The authors discuss the most critical obstacles associated with these studies, and suggest possible solutions to further advance the field.

The rapid interdisciplinary advances in nanotechnology can be characterised as the beginning of a new industrial revolution, where novel devices and materials are fabricated and controlled on the atomic level. Ion- and water-selective nanopores represent an important frontier in these advances.

The selected papers in this Special Issue provide both a snapshot of the present as well as strong indications of how the subject is likely to evolve over the coming years. We may, for example, anticipate major developments in the theory at a fundamental level, based on statistical mechanics and quantum mechanics; substantial improvements in "intermediate-level" theories like PNP, modified CKE, and DFT which promise quantitative predictions of the properties of real channels; as well as much faster and more capacious MD modelling of larger ensembles of atoms on longer timescales, more accurate due to use of polarisable force fields and QM/MM, encompassing gating and permeation events at a statistically useful level. This progress is expected to lead to the first-principles design and fabrication of structures optimised for many important applications including ion pumps, energy harvesting, and field-effect ionic transistors as well as those mentioned at the beginning. Many of these will require theory and experiment on small scales where disciplinary distinctions have mostly faded away, but where physics predominates.

An additional impulse propelling these developments forward is expected due to the fusion of physics-based approaches with artificial intelligence. The latter has already been proven to be very useful for the accelerated learning of the force fields in MD, as well as for the reconstruction of the potentials of the mean force and neural-network-based discovery of partial differential equations. Remarkably, it also underlies a recent breakthrough in the solution of the protein-folding problem.

**Author Contributions:** The authors contributed equally to this paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Leverhulme Trust (UK) grant number RPG-2017-134.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We are very grateful for much help and advice from Miroslav Barabash, Bob Eisenberg, William Gibby and other authors of the papers in this Special Issue and we acknowledge enumerable valuable discussions with Aneta Stefanovska.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Reference**


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

#### **Hamid Mofidi <sup>1</sup> , Bob Eisenberg <sup>2</sup> and Weishi Liu 1,\***


Received: 3 January 2020; Accepted: 10 March 2020; Published: 12 March 2020

**Abstract:** 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.

**Keywords:** reversal potential; effects of diffusion coefficients; permanent charge

#### **1. 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–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–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–32].

There have been great achievements in analyzing the PNP models for ionic flows through ion channels [5,28,33–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® (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:


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:


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.

#### *1.1. 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–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*,

$$\begin{aligned} \text{Poisson:} \quad & \nabla \cdot \left( \varepsilon\_{r}(\overrightarrow{\boldsymbol{X}}^{\dagger}) \varepsilon\_{0} \nabla \Phi \right) = -e\_{0} \Big( \sum\_{s=1}^{n} z\_{s} \mathbb{C}\_{s} + \mathcal{Q}(\overrightarrow{\boldsymbol{X}}^{\dagger}) \Big), \\ \text{Nernst-Planck:} \quad & \partial\_{t} \mathbb{C}\_{k} + \nabla \cdot \overrightarrow{\mathcal{J}}\_{k} = 0, \quad -\overrightarrow{\mathcal{J}}\_{k} = \frac{1}{k\_{B}T} \mathcal{D}\_{k}(\overrightarrow{\boldsymbol{X}}^{\dagger}) \mathbb{C}\_{k} \nabla \mu\_{k}. \end{aligned} \tag{1}$$

where −→*<sup>X</sup>* <sup>∈</sup> <sup>Ω</sup> with <sup>Ω</sup> being a three-dimensional cylindrical-like domain representing the channel of length *<sup>L</sup>*<sup>ˆ</sup> (nm <sup>=</sup> *<sup>L</sup>*<sup>ˆ</sup> <sup>×</sup> <sup>10</sup>−9m), <sup>Q</sup>( −→*X* ) is the permanent charge density of the channel (with unit 1M = 1Molar = 1mol/L = 103mol/m<sup>3</sup> ), *εr*( −→*X* ) is the relative dielectric coefficient (with unit 1), *<sup>ε</sup>*<sup>0</sup> <sup>≈</sup> 8.854 <sup>×</sup> <sup>10</sup>−<sup>12</sup> F m−<sup>1</sup> is the vacuum permittivity, *<sup>e</sup>*<sup>0</sup> <sup>≈</sup> 1.602 <sup>×</sup> <sup>10</sup>−19<sup>C</sup> (coulomb) is the elementary charge, *<sup>k</sup><sup>B</sup>* <sup>≈</sup> 1.381 <sup>×</sup> <sup>10</sup>−23JK−<sup>1</sup> 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−<sup>1</sup> ), and, for the *k*-th ion species, *C<sup>k</sup>* is the concentration (with unit M), *z<sup>k</sup>* is the valence (the number of charges per particle with unit 1), and *µ<sup>k</sup>* is the electrochemical potential (with unit J = CV) depending on electrical potential Φ and concentrations *C<sup>k</sup>* . The flux density −→J *<sup>k</sup>* ( −→*X* ) (with unit mol m−<sup>2</sup> s −1 ) is the number of particles across each cross-section in per unit time, D*<sup>k</sup>* ( −→*X* ) is the diffusion coefficient (with unit m2/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*,

$$\begin{split} \frac{1}{\mathcal{A}(X)} \frac{d}{dX} \left( \varepsilon\_r(X) \varepsilon\_0 \mathcal{A}(X) \frac{d\Phi}{dX} \right) &= -\varepsilon\_0 \left( \sum\_{s=1}^n z\_s \mathbb{C}\_s + \mathcal{Q}(X) \right), \\ \frac{d\mathcal{J}\_k}{dX} &= 0, \quad -\mathcal{J}\_k = \frac{1}{k\_B T} \mathcal{D}\_k(X) \mathcal{A}(X) \mathbb{C}\_k \frac{d\mu\_k}{dX}, \end{split} \tag{2}$$

where *X* is the coordinate along the channel, A(*X*) is the area of cross-section of the channel over location *X*, and J*<sup>k</sup>* (with unit mol s−<sup>1</sup> ) is the total flux through the cross-section. Equipped with System (2), we impose the following boundary conditions, for *k* = 1, 2, · · · , *n*,

$$\Phi(0) = \mathcal{V}, \quad \mathbb{C}\_{k}(0) = L\_{k} > 0; \quad \Phi(\hat{L}) = 0, \quad \mathbb{C}\_{k}(\hat{L}) = R\_{k} > 0. \tag{3}$$

One often uses the electroneutrality conditions on the boundary concentrations because the solutions are made from electroneutral solid salts,

$$\sum\_{s=1}^{n} z\_s L\_s = \sum\_{s=1}^{n} z\_s R\_s = 0.\tag{4}$$

The electrochemical potential *µ<sup>k</sup>* (*X*) for the *k*-th ion species consists of the ideal component *µ id k* (*X*) and the excess component *µ ex k* (*X*), i.e., *µ<sup>k</sup>* (*X*) = *µ id k* (*X*) + *µ ex k* (*X*). 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,

$$
\mu\_k(X) = \mu\_k^{\rm id}(X) = z\_k e\_0 \Phi(X) + k\_B T \ln \frac{\mathbb{C}\_k(X)}{\mathbb{C}\_0},\tag{5}
$$

where *C*<sup>0</sup> is a characteristic concentration of the problems, and one may consider

$$\mathcal{C}\_0 = \max\_{1 \le k \le n} \{ L\_{k'} R\_{k'} \sup\_{X \in [0, \hat{L}]} |\mathcal{Q}(X)|\}. \tag{6}$$

For given V, Q(*X*), *L<sup>k</sup>* 's and *R<sup>k</sup>* 's, if (Φ(*X*), *C<sup>k</sup>* (*X*), J*<sup>k</sup>* ) is a solution of the boundary value problem (BVP) (2) and (3), then the electric current I is

$$\mathcal{Z} = e\_0 \sum\_{s=1}^n z\_s \mathcal{J}\_s. \tag{7}$$

For an analysis of the BVP (2) and (3), we work on a dimensionless form. Set

$$\mathcal{D}\_0 = \max\_{1 \le k \le n} \{ \sup\_{X \in [0, \hat{L}]} \mathcal{D}\_k(X) \} \text{ and } \bar{\varepsilon}\_r = \sup\_{X \in [0, \hat{L}]} \varepsilon\_r(X).$$

Let

$$\begin{aligned} \varepsilon^{2} &= \frac{\mathbb{E}\_{\mathbf{r}} \varepsilon\_{0} k\_{B} T}{\varepsilon\_{0}^{2} \mathbb{L} \mathbb{C}\_{0}}, \quad \mathbb{E}\_{\mathbf{r}}(\mathbf{x}) = \frac{\varepsilon\_{\mathbf{r}}(\mathbf{X})}{\mathbb{E}\_{\mathbf{r}}}, \quad \mathbf{x} = \frac{\mathbf{X}}{\mathbf{L}}, \quad h(\mathbf{x}) = \frac{\mathcal{A}(\mathbf{X})}{\mathcal{L}^{2}}, \quad D\_{k}(\mathbf{x}) = \frac{\mathcal{D}\_{k}(\mathbf{X})}{\mathcal{D}\_{0}},\\ \mathcal{Q}(\mathbf{x}) = \frac{\mathcal{Q}(\mathbf{X})}{\mathbb{C}\_{0}}, \quad \phi(\mathbf{x}) = \frac{e\_{0}}{k\_{B}T} \Phi(\mathbf{X}), \quad c\_{k}(\mathbf{x}) = \frac{\mathcal{C}\_{k}(\mathbf{X})}{\mathcal{C}\_{0}}, \quad \mathfrak{H}\_{k} = \frac{1}{k\_{B}T} \mu\_{k}, \quad \mathfrak{J}\_{k} = \frac{\mathcal{J}\_{k}}{\mathbb{L} \mathbb{C}\_{0} \mathcal{D}\_{0}}. \end{aligned} \tag{8}$$

In terms of the new variables, the BVP (2) and (3) becomes, for *k* = 1, 2, · · · , *n*,

$$\begin{aligned} \frac{\varepsilon^2}{h(\mathbf{x})} \frac{d}{d\mathbf{x}} \left( \hat{\varepsilon}\_r(\mathbf{x}) h(\mathbf{x}) \frac{d}{d\mathbf{x}} \phi \right) &= -\sum\_{s=1}^n z\_s c\_s - Q(\mathbf{x}),\\ \frac{dJ\_k}{d\mathbf{x}} &= 0, \quad -J\_k = h(\mathbf{x}) D\_k(\mathbf{x}) c\_k \frac{d}{d\mathbf{x}} \hat{\mu}\_{k\prime} \end{aligned} \tag{9}$$

with the boundary conditions

$$\phi(0) = V = \frac{\varepsilon\_0}{k\_B T} \mathcal{V}, \quad c\_k(0) = l\_k = \frac{L\_k}{\mathbb{C}\_0}; \quad \phi(1) = 0, \quad c\_k(1) = r\_k = \frac{R\_k}{\mathbb{C}\_0}.\tag{10}$$

**Remark 1.** *The actual dimensional forms of quantities have been used for all figures throughout the paper, that is,*

$$\begin{aligned} \mathsf{C}\_{k} &= \mathsf{C}\_{0}\mathsf{c}\_{k}\ (\mathsf{M}), & \mathsf{Q} &= \mathsf{C}\_{0}\mathsf{Q}\ (\mathsf{M}), & \mu\_{k} &= \left(\mathsf{e}\_{0}\mathsf{@} + k\_{\mathsf{B}}\mathsf{T}\ln\left(\mathsf{C}\_{k}/\mathsf{C}\_{0}\right)\right)\ (\mathsf{J}),\\ \mathsf{Q} &= \frac{k\_{\mathsf{B}}\mathsf{T}}{\mathsf{e}\_{0}}\mathsf{q}\ (\mathsf{V}), & \mathcal{J}\_{k} &= \mathsf{\hat{L}}\mathsf{C}\_{0}\mathcal{D}\_{0}\mathsf{J}\_{k}\ (\mathsf{mol}/\mathsf{s}), & \mathcal{J} &= \mathsf{\hat{L}}\mathsf{C}\_{0}\mathcal{D}\_{0}\mathsf{e}\_{0}\mathsf{I}\ (\mathsf{A}), \end{aligned} \tag{11}$$

*and we take C*<sup>0</sup> <sup>=</sup> <sup>10</sup> *M, <sup>L</sup>*<sup>ˆ</sup> <sup>=</sup> 2.5 *nm and* <sup>D</sup><sup>0</sup> <sup>=</sup> 2.032 <sup>×</sup> <sup>10</sup>−<sup>9</sup> *<sup>m</sup>*2/*s*, *and, for diffusion constants [31],*

$$\begin{aligned} \mathcal{D}\_k &= 1.334 \times 10^{-9} \,\mathrm{m}^2/\mathrm{s} \,\mathrm{for}\,\mathrm{Na}^+, \,\mathrm{or} \\ \mathcal{D}\_k &= 2.032 \times 10^{-9} \,\mathrm{m}^2/\mathrm{s} \,\mathrm{for}\,\mathrm{Cl}^-, \,\mathrm{or} \\ \mathcal{D}\_k &= 0.792 \times 10^{-9} \,\mathrm{m}^2/\mathrm{s} \,\mathrm{for}\,\,\mathrm{Ca}^{2+}. \end{aligned} \tag{12}$$

#### *1.2. 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


$$Q(\mathbf{x}) = \begin{cases} Q\_1 = Q\_3 = 0, & \mathbf{x} \in (0, a) \cup (b, 1), \\ Q\_2 = 2Q\_{0\prime} & \mathbf{x} \in (a, b), \end{cases} \tag{13}$$

where *Q*<sup>0</sup> is a constant.

We assume that *ε* > 0 in System (14) is small. The assumption is reasonable since, if *L*ˆ = 2.5 nm = 2.5 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup> and *<sup>C</sup>*<sup>0</sup> <sup>=</sup> <sup>10</sup> M, then *<sup>ε</sup>* <sup>≈</sup> <sup>10</sup>−<sup>3</sup> [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,

$$\begin{aligned} \varepsilon \phi &= u, \quad \varepsilon \dot{u} = -\sum\_{s=1}^{2} z\_s c\_s - Q(\mathbf{x}) - \varepsilon \frac{h\_k(\mathbf{x})}{h(\mathbf{x})} u, \\ \varepsilon \dot{c}\_k &= -z\_k c\_k u - \varepsilon \frac{f\_k}{D\_k h(\mathbf{x})}, \quad \dot{f}\_k = 0. \end{aligned} \tag{14}$$

System (14) will be treated as a dynamical system with the phase space R<sup>7</sup> 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*<sup>0</sup> where *l*<sup>1</sup> = *l*<sup>2</sup> = *l*, *r*<sup>1</sup> = *r*<sup>2</sup> = *r*, *Q*<sup>2</sup> = 2*Q*0.

**Remark 2.** *We remind the readers that the quantities V*, *l*,*r*, *c<sup>k</sup>* , *Q*, *φ*, *µ*ˆ *<sup>k</sup>* , *Jk* , *D<sup>k</sup> , and I are dimensionless quantities corresponding to the dimensional quantities* V, *L*, *R*, *C<sup>k</sup>* , Q, Φ, *µ<sup>k</sup>* , J*<sup>k</sup>* , D*<sup>k</sup> , 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.*

#### **2. 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<sup>1</sup> = D<sup>2</sup> and where D<sup>1</sup> 6= D2, 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 *∂tc* = D*∂* 2 *xxc*, 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–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*<sup>1</sup> = −*z*2. It follows that the reduced matching system for zero current *I* = *J*<sup>1</sup> − *J*<sup>2</sup> = 0 when *z*<sup>1</sup> = −*z*<sup>2</sup> = 1 is

$$\mathcal{G}\_1(A, \mathcal{Q}\_{0\prime}\theta) = V \text{ and } \mathcal{G}\_2(A, \mathcal{Q}\_{0\prime}\theta) = 0,\tag{15}$$

where

$$\begin{aligned} G\_1(A, Q\_0, \theta) &= \theta \left( \ln \frac{S\_a + \theta Q\_0}{S\_b + \theta Q\_0} + \ln \frac{l}{r} \right) - (1 + \theta) \ln \frac{A}{B} + \ln \frac{S\_a - Q\_0}{S\_b - Q\_0}, \\ G\_2(A, Q\_0, \theta) &= \theta Q\_0 \ln \frac{S\_a + \theta Q\_0}{S\_b + \theta Q\_0} - N\_\prime \end{aligned} \tag{16}$$

and, *A* is the geometric mean of concentrations at *x* = *a*, that is,

$$A = \sqrt{c\_1(a)c\_2(a)}\tag{17}$$

$$B = \frac{1 - \beta}{a}(l - A) + r, \quad S\_d = \sqrt{Q\_0^2 + A^2}, \quad S\_b = \sqrt{Q\_0^2 + B^2}, \quad N = A - l + S\_d - S\_{b\nu} \tag{18}$$

and

$$\theta = \frac{D\_2 - D\_1}{D\_2 + D\_1}, \quad a = \frac{H(a)}{H(1)}, \quad \beta = \frac{H(b)}{H(1)}\text{ where }H(\mathbf{x}) = \int\_0^\mathbf{x} \frac{1}{h(\mathbf{s})} d\mathbf{s}.\tag{19}$$

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*<sup>0</sup> 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.

For fixed *Q*<sup>0</sup> and *θ*, *A* can actually be solved from *G*2(*A*, *Q*0, *θ*) = 0, where *G*<sup>2</sup> is defined in Display (16) with the properties stated in the next theorem.

**Theorem 1** (Theorem 3.4 in [23])**.** *The solution A* = *A*(*Q*0, *θ*) *of G*2(*A*, *Q*0, *θ*) = 0 *satisfies*


*(d) if θQ*<sup>0</sup> ≥ 0*, then ∂Q*<sup>0</sup> *A*(*Q*0, *θ*) *has the same sign as that of* (*l* − *r*)*Q*0*,*

*where A*∗ = (1 − *β*)*l* + *αr* 1 − *β* + *α .*

For fixed *Q*<sup>0</sup> and *θ*, the reversal potential *Vrev* = *Vrev*(*Q*0, *θ*) can also be determined and enjoy properties stated in the next two theorems. Recall that we denote *J*<sup>1</sup> = *J*<sup>2</sup> by *J*.

**Theorem 2** (Theorem 4.2 in [23])**.** *For the reversal potential Vrev* = *Vrev*(*Q*0, *θ*)*, one has*

$$\text{(i)}\quad \text{if } l > r \text{, then } l > 0 \text{, and, hence, } -\frac{1}{z\_1}\ln\frac{l}{r} < V\_{\text{rev}}(Q\_0, \theta) < \frac{1}{z\_1}\ln\frac{l}{r};$$

$$\text{If (ii)} \quad \text{if } l < r \text{, then } \mathbf{J} < 0 \text{, and, hence, } \frac{1}{z\_1} \ln \frac{l}{r} < V\_{rev}(Q\_{0r}\theta) < -\frac{1}{z\_1} \ln \frac{l}{r}$$

$$(\text{iii)}\quad V\_{\text{rev}}(0,\theta) = \frac{\theta}{z\_1}\ln\frac{l}{r}\,and\,\lim\_{Q\_0\to\pm\infty}V\_{\text{rev}}(Q\_0,\theta) = \pm\frac{1}{z\_1}\ln\frac{l}{r}.$$

**Theorem 3** (Theorem 4.4 in [23])**.** *For any given θ* ∈ (−1, 1)*, one has*


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.

#### *2.1. Zero-Current Flux J* = *J*<sup>1</sup> = *J*2*.*

We aim to clarify the relationships of ion fluxes with permanent charge and diffusion constants when current is zero.

Recall that fluxes *J*<sup>1</sup> and *J*<sup>2</sup> are equal for this case and let *J* denote it. For any permanent charge *Q* = 2*Q*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

$$J = -\frac{6D\_1 D\_2 (A - l)}{(D\_1 + D\_2)} = -\frac{6D\_1 D\_2 (r - B)}{(D\_1 + D\_2)}.\tag{20}$$

*;*

## 2.1.1. 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,

$$\frac{J\_k}{D\_k} \int\_0^1 \frac{1}{h(\mathbf{x})c\_k(\mathbf{x})} d\mathbf{x} = z\_k V + \ln \frac{l}{r}.\tag{21}$$

Therefore, the sign of flux *J<sup>k</sup>* 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* = *Vrev* 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 *<sup>l</sup>* <sup>−</sup> *<sup>r</sup>*. This is consistent with observations in Figure <sup>1</sup> where <sup>D</sup><sup>1</sup> <sup>=</sup> 1.334 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>2/s is fixed, and <sup>D</sup><sup>2</sup> varies from the same value to <sup>D</sup><sup>2</sup> <sup>=</sup> 2.032 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>2/s, and to a random large value.

**Figure 1.** The function J = J (Q) for various values of *ρ* = *D*2/*D*<sup>1</sup> : The left panel for *L* = 2 mM and *R* = 5 mM; the right panel for *L* = 5 mM and *R* = 2 mM.

2.1.2. Dependence of Zero-Current Flux *J* on *Q*<sup>0</sup> and *D<sup>k</sup>* 's

Concerning the dependence of the zero-current flux *J* on *Q*0, we have the following:

(i) *If D*<sup>1</sup> = *D*2*, then the zero-current flux J is an even function in Q*0*, and it is monotonic for Q*<sup>0</sup> > 0*.*

In this case, *θ* = 0 and, hence, it follows from Theorem 1 that *A* is an even function in *Q*<sup>0</sup> and is monotonic in *Q*<sup>0</sup> for *Q*<sup>0</sup> > 0, and thus is the zero-current flux *J* from Display (20).

(ii) *If D*<sup>1</sup> 6= *D*2*, then the zero-current flux J is not an even function in Q*<sup>0</sup> *and the monotonicity of the zero-current flux J in Q*<sup>0</sup> *seems to be more complicated.*

In this case, it can be seen that *G*<sup>2</sup> 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*<sup>0</sup> → ±∞ that is consistent with the observations in Figure 1.

(iii) Another fascinating result is that *the magnitude of ρ* = *D*2/*D*<sup>1</sup> *affects the monotonicity of the zero-current flux J in Q*0*.*

In this case, if one fixes *D*1, and let *D*<sup>2</sup> increases from small values to *D*<sup>2</sup> → ∞, (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*<sup>0</sup> that is not intuitive.

Let us consider the case where *L* < *R* and *Q*<sup>0</sup> < 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*<sup>0</sup> increases,


Thus, depending on the size of *ρ*, the zero-current flux *J* may increase or decrease in *Q*<sup>0</sup> < 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.

#### *2.2. Reversal Potential Vrev.*

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., <sup>V</sup> <sup>=</sup> <sup>Φ</sup>(0) <sup>−</sup> <sup>Φ</sup>(*L*ˆ) that produces zero current <sup>I</sup>. Using dimensionless form of quantities (see Remark 2), it follows from System (15) and (16) (where there are two ion species with valences *z*<sup>1</sup> = −*z*<sup>2</sup> = 1) that for general permanent charge *Q* = 2*Q*<sup>0</sup> 6= 0 with arbitrary diffusion constants [23], the variable *A* (the geometric mean of concentrations at *x* = *a*) can be solved uniquely from *G*<sup>2</sup> = 0 in System (15), and the reversal potential is then

$$V\_{rev} = \theta \left( \ln \frac{S\_d + \theta Q\_0}{S\_b + \theta Q\_0} + \ln \frac{l}{r} \right) - (1 + \theta) \ln \frac{A(Q\_0, \theta)}{B(Q\_0, \theta)} + \ln \frac{S\_d - Q\_0}{S\_b - Q\_0} \tag{22}$$

where *B*, *Sa*, *S<sup>b</sup>* , and *θ* are defined in Displays (18) and (19).

#### 2.2.1. Range of Reversal potential *Vrev*

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

#### 2.2.2. Zero Reversal Potential

One particular case is when the reversal potential is zero. To examine under what conditions one obtains *Vrev* = *Vrev*(*Q*0) = 0, it follows Theorem 2 that,


Considering the second case above, the observations in Figure 2 show that, as *ρ* = *D*2/*D*<sup>1</sup> increases, magnitude of the corresponding *Q*<sup>0</sup> becomes larger. In fact, as *ρ* → ∞, then *Q*<sup>0</sup> → −∞.

**Figure 2.** The function V = V*rev*(Q): The left panel for *L* = 2 mM and *R* = 5 mM; the right panel for *L* = 5 mM and *R* = 2 mM.

#### 2.2.3. Reversal Potential *Vrev*(*Q*0) for *Q*<sup>0</sup> = 0

For *Q*<sup>0</sup> = 0, one has *Vrev*(0) = *θ* ln *<sup>l</sup> r* from Theorem 2. Therefore,


Let us consider the case where *D*<sup>1</sup> < *D*2. In that case, *Vrev*(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 *Jk*/*D<sup>k</sup>* are independent of *D<sup>k</sup>* in Equation (21)), and to help |*J*1| more than |*J*2| to get *J*<sup>1</sup> = *J*<sup>2</sup> for zero current conditions, one needs to increase *V* when *l* > *r* and decrease *V* when *l* < *r* , and that is why *Vrev*(0) > 0 for *l* > *r* and *Vrev*(0) < 0 for *l* < *r* . This is consistent with observations in Figure 2 as well. The analysis for the other case with *D*<sup>1</sup> > *D*<sup>2</sup> is similar.

2.2.4. Monotonicity of V*rev* with Respect to Q

It follows from Theorem 3 that

*For θQ* > 0, *∂*QV*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* : *Vrev is increasing in Q f or l* > *r and decreasing in Q f or l* < *r*.

We remark that, in Figure 3, we take *<sup>L</sup>* <sup>=</sup> <sup>20</sup> mM, *<sup>R</sup>* <sup>=</sup> <sup>50</sup> mM, and <sup>D</sup><sup>1</sup> <sup>=</sup> 1.334 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>2/s and <sup>D</sup><sup>2</sup> <sup>=</sup> 2.032 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>2/s which are diffusion constants of, say, Na<sup>+</sup> and Cl−, respectively (see the solid line), and <sup>D</sup><sup>1</sup> <sup>=</sup> 1.334 <sup>×</sup> <sup>10</sup>−9m2/s and <sup>D</sup><sup>2</sup> <sup>=</sup> 0.792 <sup>×</sup> <sup>10</sup>−9m2/s, where <sup>D</sup><sup>2</sup> is the diffusion constants of Ca2<sup>+</sup> (see the dashed line).

**Figure 3.** V = V*rev*(Q) decreases when *L* < *R*, independent of values of diffusion constants.

2.2.5. Dependence of V*rev* on *ρ* = *D*2/*D*<sup>1</sup>

Let us discuss the dependence of *Vrev* on *ρ* = *D*2/*D*<sup>1</sup> for effects of *D*<sup>1</sup> and *D*2. It follows from Proposition 4.6 in [23] that

#### *The reversal potential Vrev 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<sup>k</sup>* are independent of *D<sup>k</sup>* [36]. The relation surely holds for the zero-current condition, i.e., *J*<sup>1</sup> = *J*<sup>2</sup> with *V* = *Vrev*. Now, let us fix *D*<sup>1</sup> and increase *D*<sup>2</sup> (so *ρ* is increasing). Then, |*J*2| increases since all but *J*2/*D*<sup>2</sup> in Equation (21) are independent of *D*2. Consequently, to meet the zero-current condition, we need to increase |*J*1|. Intuitively increasing *Vrev* 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, *Vrev* decreases. This counterintuitive behavior

could be explained by the fact that *c*1(*x*) depends on *Vrev*, and reducing *Vrev* could still increase |*J*1|. In fact, *l* < *r* will result in reducing *Vrev*, 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, *<sup>L</sup>*<sup>ˆ</sup> and <sup>D</sup><sup>1</sup> for Na<sup>+</sup> as mentioned in Remark 1. Now, suppose that <sup>D</sup><sup>1</sup> <sup>2</sup> <sup>=</sup> 0.792 <sup>×</sup> <sup>10</sup>−9m2/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 <sup>J</sup> <sup>=</sup> <sup>−</sup>1.7632 <sup>×</sup> <sup>10</sup>−17mol s−<sup>1</sup> . Now, if we increase <sup>D</sup><sup>1</sup> 2 to <sup>D</sup><sup>2</sup> <sup>2</sup> <sup>=</sup> 2.032 <sup>×</sup> <sup>10</sup>−9m2/s, which is Cl<sup>−</sup> diffusion constant, then <sup>V</sup>*rev* <sup>=</sup> <sup>−</sup>19.5527 mV and <sup>J</sup> <sup>=</sup> <sup>−</sup>1.8788 <sup>×</sup> <sup>10</sup>−17mol s−<sup>1</sup> . 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<sup>1</sup> is fixed, but we increase D2.

## *2.3. 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

$$I = V \sum\_{k=1}^{n} z\_k^2 D\_k \frac{r\_k - l\_k e^{z\_k V}}{1 - e^{z\_k V}}.$$

For the case where *n* = 2 and *z*<sup>1</sup> = −*z*<sup>2</sup> = 1, the GHK equation for the reversal potential is

$$V\_{rev}^{GHK}(\rho) = \ln \frac{r + \rho l}{l + \rho r}. \tag{23}$$

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*<sup>0</sup> = 0. However, this is not correct either. It follows from Formula (22) for *Q*<sup>0</sup> = 0 that the zeroth order in *ε* approximation of the reversal potential in this case is

$$V\_{rev}(0,\rho) = \frac{\rho - 1}{\rho + 1} \ln \frac{l}{r}.\tag{24}$$

Figure 5 compares *Vrev*(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.

**Figure 5.** V*rev*(*Q* = 0, *ρ*) vs. V *GHK rev* (*ρ*): The left panel for *L* = 20 mM and *R* = 50 mM; the right panel for *L* = 1 mM and *R* = 50 mM.

In Figure 6, we arrange a simple numerical result for the case where Q 6= 0 to compare the graphs of V*rev*(Q, *ρ*), obtained from Formula (22), for various values of permanent charge Q. We consider *L* = 20 mM, *R* = 50 mM, and 0 < *ρ* < 5 for some values of Q, i.e., Q = 0 M, 1 M, 10 M.

**Figure 6.** V*rev*(Q, *ρ*) with various values of permanent charges.

#### *2.4. 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*, D1, D2) to numerically investigate the behavior of *C<sup>k</sup>* (*X*) and Φ(*X*) throughout the channel. Figures 7 and 8 graph the cases with small permanent charge Q = 0.1 mM when *L* = 20 mM, *R* = 50 mM, D<sup>1</sup> =

1.334 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>2/s, and <sup>D</sup><sup>2</sup> <sup>=</sup> 2.032 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>2/s. In this case, we obtain <sup>J</sup> <sup>=</sup> <sup>−</sup>1.2079 <sup>×</sup> <sup>10</sup>−<sup>16</sup> mol s−<sup>1</sup> and V*rev* = −4.4820 mV.

**Figure 7.** The functions *C<sup>k</sup>* (*X*) (**left**) and Φ(*X*) (**right**) with Q = 0.1 mM.

**Figure 8.** The functions *µ*1(*X*) and *µ*2(*X*) are increasing for Q = 0.1 mM.

Furthermore, Figures 9 and 10 show graphs of concentrations, electrical potential, and electrochemical potentials versus *X*, where *L* = 20 mM, *R* = 50 mM, *Q* = 2 M, and diffusion constants are the same as the previous one. In this case, we obtain <sup>J</sup> <sup>=</sup> <sup>−</sup>1.8789 <sup>×</sup> <sup>10</sup>−<sup>17</sup> mol s−<sup>1</sup> and V*rev* = −19.5527 mV.

**Remark 3.** *We end this section with a few of the remarks on some important features captured in Figures 7–10. It follows from the Nernst–Planck equation that µ* 0 *k* (*x*) *has the same sign as that of µ<sup>k</sup>* (1) − *µ<sup>k</sup>* (0) *or the opposite sign as that of J<sup>k</sup> ; in particular, µ<sup>k</sup>* (*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 µ<sup>k</sup>* (*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, µ* 0 *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, µ* 0 *k* (*x*) > 0 *for k* = 1, 2*, and, hence, µ<sup>k</sup>* (*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<sup>k</sup>* (*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.*

**Figure 9.** The functions *C*1(*X*) and *C*2(*X*) (**left**) and the function Φ(*X*) (**right**) for Q = 1 M.

**Figure 10.** The functions *µ*1(*X*) and *µ*2(*X*) are increasing for Q = 1 M.

#### **3. 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<sup>k</sup>* (*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.

#### *3.1. 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

$$
\eta \ln \frac{S\_b - \eta}{S\_a - \eta} - N = 0,\tag{25}
$$

where *B* = *l* − *A* + *r*, *Sa*, *S<sup>b</sup>* and *N* are defined in Display (18), and

$$\eta = Q\_0 - \frac{Q\_0}{\ln \frac{Bl}{Ar}} \left( V + \ln \frac{l(S\_b - Q\_0)}{r(S\_a - Q\_0)} \right) + \frac{N}{\ln \frac{Bl}{Ar}}.\tag{26}$$

Once *A* is solved from Equation (25), we can obtain the flux densities and current equations as follows:

$$\begin{aligned} \mathcal{J}\_k &:= \mathcal{J}\_k(V, l, r, D\_1, D\_2) = 3D\_k(l-A) \left( 1 + (-1)^k \frac{\eta}{Q\_0} \right), \quad \text{for } k = 1, 2, \\\ I &:= I(V, l, r, D\_1, D\_2) = \mathcal{J}\_1 - \mathcal{J}\_2 = 3(l-A) \left( D\_1 - D\_2 - \frac{\eta}{Q\_0} (D\_1 + D\_2) \right). \end{aligned} \tag{27}$$

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.

#### *3.2. Current–Voltage and Current-Permanent Charge Relations*

#### 3.2.1. 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) that, if *<sup>V</sup>* <sup>=</sup> ln *<sup>B</sup>*(*S<sup>a</sup>* <sup>−</sup> *<sup>Q</sup>*0) *A*(*S<sup>b</sup>* − *Q*0) , then *η* = 0. Therefore, from Display (27), *I* = 3(*l* − *A*)(*D*<sup>1</sup> − *D*2), which implies

*For special values o f parameters* (*l*,*r*, *V*, *Q*0), *the sign o f I is determined by the sign o f D*<sup>1</sup> − *D*2.

#### 3.2.2. I–V Curves and I–Q Curves

Figure 11 is a numerical simulation from Equations (25) and (26) of the I–V curves for several values of <sup>Q</sup> with <sup>D</sup><sup>1</sup> <sup>=</sup> 1.334 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>2/s and <sup>D</sup><sup>2</sup> <sup>=</sup> 2.032 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup>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].

**Figure 11.** The function I = I(V) for *L* = 20 mM and *R* = 50 mM.

Now, for I–Q relations, our numerical experiments shows that: *I–Q curves are not monotonic in general.* Recall that Equation (21), in dimensional form, gives

$$\mathcal{I}\_k \int\_0^\mathbb{L} \frac{k\_B T}{\mathcal{D}\_k \mathcal{A}(X) \mathbb{C}\_k(X)} dX = \mu\_k(0) - \mu\_k(\hat{L}), \quad k = 1, 2.$$

The sign of J*<sup>k</sup>* is determined by the boundary conditions, independent of the permanent charge. Nevertheless, as expected and seen in Figure 12, the magnitudes of J*<sup>k</sup>* 's, and, consequently, the sign and the size of the current I do depend on Q = 2*Q*0*C*<sup>0</sup> 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,


In particular, I–Q curves are not monotonic in general.

**Figure 12.** The function I = I(Q) with D<sup>1</sup> < D2: The left panel for *L* = 20 mM and *R* = 30 mM; the right panel for *L* = 30 mM and *R* = 20 mM.

In addition, we claim based on numerical observations (not proven though) that there exists <sup>V</sup>ˆ(D1, <sup>D</sup>2) = min{V<sup>∗</sup> , |V∗|}, such that


In particular, it can be seen in Section 3.3 that current is monotonic in Q for V = 0. In the end, we would like to mention that the diffusion constants affect the values *V* <sup>∗</sup> and *V*<sup>∗</sup> above.

#### *3.3. Zero-Voltage Current*

The different permeability of the membrane determines the zero membrane potential (voltage) to different types of ions, as well as the concentrations of the ions, the permanent charge, and the shape of the channel. Denote current *I*(*V*; *Q*), and the fluxes *J<sup>k</sup>* (*V*; *Q*), for *k* = 1, 2, to include the dependence on *Q* too. We call *I*(0; *Q*) the *zero-potential current* and *J<sup>k</sup>* (0; *Q*) the *zero-potential fluxes*, respectively, when *V* = 0. For any given value of membrane potential *V*, approximation formulas for the current *I*(*V*; *Q*), for small and large values of permanent charge *Q*, are provided in [34,38], respectively.

It follows from [34] that, for *small* values of *Q*, applying *V* = 0, zero-potential current *I s* (0; *Q*), and zero-potential fluxes *J s k* (0; *Q*) (in dimensionless forms as mentioned in Remark 2) are

$$\begin{split} I^s(0;Q) &= (l-r)\left(D\_1 - D\_2\right) - \frac{3}{2} \frac{(D\_1 + D\_2)(l-r)^2}{(2l+r)(l+2r)\ln\frac{l}{r}} Q + O(Q^2), \\ I^s\_k(0;Q) &= (l-r)D\_k + (-1)^k \frac{3(l-r)^2 D\_k}{2(2l+r)(l+2r)\ln\frac{l}{r}} Q + O(Q^2), \quad k = 1,2. \end{split} \tag{28}$$

For *large* positive values of *Q* = 2*Q*0, with *ν* = 1/*Q*<sup>0</sup> (where *ν* is small), it follows from [38] that zero-potential current *I l* (*ν*) = *I l* (0; *Q*) and zero-potential fluxes *J l k* (*ν*) = *J l k* (0; *Q*) are

$$\begin{split} I^{l}(\nu) &= -6D\_{2}\sqrt{lr}\frac{(\sqrt{l}-\sqrt{r})}{\sqrt{l}+\sqrt{r}} + \frac{3}{2}D\_{1}\left(\frac{l+r}{\sqrt{l}+\sqrt{r}}\right)^{2}(l-r)\nu \\ &+ \frac{3}{2}D\_{2}\frac{l+r}{(\sqrt{l}+\sqrt{r})^{2}}f(l,r)(l-r)\nu + O(\nu^{2}), \\ J\_{1}^{l}(\nu) &= \frac{3}{2}D\_{1}\left(\frac{l+r}{\sqrt{l}+\sqrt{r}}\right)^{2}(l-r)\nu, \\ J\_{2}^{l}(\nu) &= 6D\_{2}\sqrt{lr}\frac{\sqrt{l}-\sqrt{r}}{\sqrt{l}+\sqrt{r}} - \frac{3}{2}D\_{2}\frac{l+r}{(\sqrt{l}+\sqrt{r})^{2}}f(l,r)(l-r)\nu + O(\nu^{2}), \end{split} \tag{29}$$

where

$$f(l,r) = \frac{2lr}{(\sqrt{l} + \sqrt{r})^2} + l + r - \frac{1}{2} \frac{\ln l - \ln r}{l - r} (l + r)^2. \tag{30}$$

It can be readily seen from Equations (28) that, for small values of *Q*, the zero-potential current *I s* (0; *Q*) is increasing in *Q* when *l* < *r* and is decreasing in *Q* if *l* > *r*.

However, for large values of the permanent charge *Q*, the zero-potential current *I l* (0; *Q*) depends on *Q* in a much richer way. To state the results, we need the following lemma.

**Lemma 1.** *There are t*<sup>1</sup> *and t*<sup>2</sup> *with* 0 < *t*<sup>1</sup> < 1 < *t*<sup>2</sup> *so that f*(*l*,*r*) > 0 *for l*/*r* ∈ (*t*1, *t*2) *and f*(*l*,*r*) < 0 *for l*/*r* 6∈ [*t*1, *t*2]*.*

Note that

$$\frac{d}{d\nu}I^l(0) = \frac{\Im D\_2}{2} \left(\frac{l+r}{\sqrt{l}+\sqrt{r}}\right)^2 \left(\frac{D\_1}{D\_2} + \frac{f(l,r)}{l+r}\right)(l-r).$$

It follows from Equations (29) and Lemma 1 that, for large values of *Q* (small values of *ν*),


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* 6= 0, then the current *I* is not monotonic in *Q*.

**Figure 13.** The function I = I(Q) for V = 0: The left panel for *L* = 20 mM and *R* = 30 mM; the right panel for *L* = 30 mM and *R* = 20 mM.

#### **4. 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.

**Author Contributions:** All authors of this paper contributed equally. All authors have read and agreed to the published version of the manuscript.

**Acknowledgments:** We thank the anonymous referees for suggestions that significantly improved the paper. The research is partially supported by Simons Foundation Mathematics and Physical Sciences-Collaboration Grants for Mathematicians #581822.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Review*
