Molecular Mean-Field Theory of Ionic Solutions: A Poisson-Nernst-Planck-Bikerman Model

We have developed a molecular mean-field theory—fourth-order Poisson–Nernst–Planck–Bikerman theory—for modeling ionic and water flows in biological ion channels by treating ions and water molecules of any volume and shape with interstitial voids, polarization of water, and ion-ion and ion-water correlations. The theory can also be used to study thermodynamic and electrokinetic properties of electrolyte solutions in batteries, fuel cells, nanopores, porous media including cement, geothermal brines, the oceanic system, etc. The theory can compute electric and steric energies from all atoms in a protein and all ions and water molecules in a channel pore while keeping electrolyte solutions in the extra- and intracellular baths as a continuum dielectric medium with complex properties that mimic experimental data. The theory has been verified with experiments and molecular dynamics data from the gramicidin A channel, L-type calcium channel, potassium channel, and sodium/calcium exchanger with real structures from the Protein Data Bank. It was also verified with the experimental or Monte Carlo data of electric double-layer differential capacitance and ion activities in aqueous electrolyte solutions. We give an in-depth review of the literature about the most novel properties of the theory, namely Fermi distributions of water and ions as classical particles with excluded volumes and dynamic correlations that depend on salt concentration, composition, temperature, pressure, far-field boundary conditions etc. in a complex and complicated way as reported in a wide range of experiments. The dynamic correlations are self-consistent output functions from a fourth-order differential operator that describes ion-ion and ion-water correlations, the dielectric response (permittivity) of ionic solutions, and the polarization of water molecules with a single correlation length parameter.

There is another important property that classical continuum theories fail to describe, namely short-range ion-ion or ion-water correlations in ion channels [8,9], charge-induced thickening and density oscillations near highly charged surfaces [14], correlation-induced charge inversion on macroions (DNA, actin, lipid membranes, colloidal particles) [59], the phase structure of plasma and polar fluids [60], colloidal charge renormalization [60], etc. Several other properties related to correlations such as the dielectric response of electrolytes solutions and the polarization of water in various conditions or external fields are usually modeled differently from the correlation perspective [61][62][63].
We have recently developed a molecular mean-field theory called-Poisson-Nernst-Planck-Bikerman (PNPB) theory-that can describe the size, correlation, dielectric, and polarization effects of ions and water in aqueous electrolytes at equilibrium or nonequilibrium all within a unified framework [64][65][66][67][68][69][70][71][72][73][74][75][76]. Water and ions in this theory can have different shapes and volumes necessarily with intermolecular voids. The theory generalizes and unifies the second-order Poisson-Bikerman equation [42] of binary ionic liquids for different-sized ions with identical steric energies [72] and the fourth-order differential permittivity operator in Santangelo's model of one component plasma [77] or in the Bazant, Storey, and Kornyshev theory of general nonlocal permittivity for equal-sized ions in ionic liquids [78].
Ion-ion and ion-water correlations are modeled by the permittivity operator with a correlation length that depends on the diameter of ions or water and the valence of ions of interest [78]. The fourth-order operator yields a permittivity as an output function of spatial variables, salt concentration, and hydration shell structure including water diameter from solving the PNPB model and thus describes the dehydration of ions from bath to channel pore or from bulk to charged wall, the polarization of water, and the change of permittivities of electrolyte solutions at different locations in response to different configurations and conditions. Water densities also change with configurations and conditions. The fourth-order operator introduces correlations into the mean-field equations so they can deal more realistically with real systems in which the correlations are of the greatest importance. A remark should be made here that simulations containing only particles do not automatically deal with correlations better than mean-field theories with fourth-order operators like this. It is not at all clear that simulations widely done in biophysics actually compute correlations well. Indeed, it is difficult to see how simulations that use conventions to approximate the electric field, and periodic boundary conditions to approximate macroscopic systems could deal with correlations correctly. The dearth of direct checks of the role of periodic boundary conditions, and of the accuracy of the conventional treatment of electrostatics, does little to assuage these concerns. The detailed direct checks found necessary in computational electronics are not easily found in simulations of ions in electrolyte solutions (see Chap. 6, particularly Figures 6.34-35 of [55] for some details that are found to be necessary in the simulations of computational electronics).
It is important to reiterate the obvious. Our model includes water as a molecule and depends on the hydration structure around ions. Our model uses partial differential equations (PDEs) to describe these essentially discrete properties of ionic solutions, and uses the physical parameters of individual atoms and water molecules, NOT just their mean-field description. This use of PDEs to describe inherently discrete processes is hardly new: most of probability theory [79,80] and the entire theory of wave equations, including the wave equation of the electron called the Schrödinger equation [81], treat discrete processes the same way, using PDEs that measure (in probability theory) the underlying discrete system, or represent it exactly as the discrete solutions of a continuum PDE (e.g., the Schrödinger equation describing a hydrogen atom).
The most important contribution of our work is to include water as discrete molecules by using Fermi distributions [82] of classical particles having excluded volumes with interstitial voids. We show that the treatment of water as finite size molecules requires, as a matter of mathematics, not physics, the existence of voids. This is demonstrated by mathematics and simple ways to compute the voids and their role are presented. These Fermi-like distributions yield saturation of all particles (ions and water) even under mathematically infinite large external fields. These distributions also satisfy mass conservation in the region of interest such as channel pores, which classical theories fail to describe as well. This Fermi distribution of classical particles obeying volume exclusion is reminiscent of the Fermi distribution of identical particles obeying the Pauli exclusion principle [83] in quantum mechanics.
We also introduce a new concept of distance-dependent potential between non-bonded particles for different-sized particles similar to the electric potential for different-charged ions and name it the steric potential. The void distribution function describes the van der Waals potential [84] of paired particles [85,86] in the system in a mean-field sense. The steric potential can be written as a distribution function of voids, emphasizing the crucial role of voids in our theory. The specific sizes of particles and the distance-dependent steric potential allow us to calculate steric energies at the atomic scale. Using Coulomb's law allows calculation of electric energies at the atomic scale as well. Therefore, our theory applies to biological or chemical systems with explicit atomic structures, as well as classical mean-field representations of bulk solutions, for example. We have shown that solving the PNPB model in different continuum and molecular domains yields self-consistent electric and steric potentials in many examples of biological ion channels or chemical systems in [64][65][66][67][68][69][70][71][72][73][74][75][76]. The theory is also consistent with classical theories in the sense that its model converges to the corresponding classical one when the volume of all particles and the correlation length tend to zero, i.e., steric and correlation effects vanish asymptotically to classical cases.
In this review article, we explain the above bold-face terms in detail and compare them with those of earlier theories in a precise but limited way. The precision means that we display explicitly, to the best of our ability, the significant differences between analogous concepts in our theory and previous treatments. It is obviously impossible to do complete comparisons in this vast and formidable field. No doubt we are ignorant of significant relevant papers. We apologize to those inadvertently slighted and ask them to help us remedy our oversight. The remaining of this article consists as follows.
Section 2 describes the physical meaning of Fermi distributions and the steric potential of ions and water with excluded volumes. We also explain the differences between Fermi and Boltzmann distributions in the context of statistical thermodynamics. Section 3 unifies Fermi distributions and correlations into the simple and concise 4th-order Poisson-Bikerman (4PBik) equation. The simplicity refers to the correlation length being the only empirical parameter in the equation. The conciseness means that the fourth-order differential operator can describe the complex and correlated properties of ion-ion and ion-water interactions, polarization of water, and dielectric response of electrolytes solutions all in a single model setting.
Section 4 presents a Gibbs free energy functional for the 4PBik equation. We show that minimization of the functional yields the equation and Fermi distributions that reduce to Boltzmann distributions when the volumes of particles vanish in limiting case. This functional is critical to explain a major shortcoming of earlier modified PB models that cannot yield Boltzmann distributions in the limit. These models are thus not consistent with classical theories and may poorly estimate steric energies and other physical properties due to their coarse approximation of size effects.
Section 5 generalizes the 4PBik equation to the PNPB model to describe flow dynamics of ions and water in the system subject to external fields. The most important feature in this section is the introduction of the steric potential to the classical Nernst-Planck equation. Electric and steric potentials describe the dynamic charge/space competition between ions and water. We also show that the PNPB model reduces to the 4PBik equation at equilibrium. Section 6 presents a generalized Debye-Hückel theory from the 4PBik equation for thermodynamic modeling. The theory yields an equation of state that analytically models ion activities in all types of binary and multi-component electrolyte solutions over wide ranges of concentration, temperature, and pressure. It is also useful to study the size, correlation, dielectric, and polarization effects in a clear comparison with those ignoring these effects. Section 7 discusses numerical methods for solving the PNPB model that is highly nonlinear and complex when coupled with the electrical field generated by protein charges in ion channels, for example. It is very challenging to numerically solve the model with tolerable accuracy in 3D protein structures that generate extremely large electric field, e.g., 0.1 V in 1 Angstrom, in parts of the molecule of great biological importance where crowded charges directly control biological function, in the same sense that a gas pedal controls the speed of a car. Section 8 demonstrates the usefulness of the PNPB theory for a wide range of biological and chemical systems, where the steric and correlation effects are of importance. We choose a few examples of these systems, namely electric double layers, ion activities, and biological ion channels.
Section 9 summarizes this review with some concluding remarks.

Fermi Distributions and Steric Potential
The total volume of an aqueous electrolyte system with K species of ions in a solvent domain Ω s is where K + 1 and K + 2 denote water and voids, respectively, v i is the volume of each species i particle, N i is the total number of species i particles, and V K+2 is the total volume of all the voids [68]. The volume of each particle v i will play a central role in our analysis, as well that the limit v i goes to zero. This limit defines the solution of point particles of classical PB and PNP theory. We must include the voids as a separate species if we treat ions and water having volumes in a model. This necessity can be proven by mathematics (see below). It is also apparent to any who try to compute a model of this type with molecular water, as it was to us [68]. Dividing the volume Equation (1) in bulk conditions by V, we get the bulk volume fraction of voids where C B i = N i V are bulk concentrations. If the system is spatially inhomogeneous with variable electric or steric fields, as in realistic systems, the constants C B i then change to functions C i (r) and so does Γ B to a void volume function We define the concentrations of particles (i.e., the distribution functions of the number density) in Ω s [72] as where φ(r) is an electric potential, S trc (r) is called a steric potential, β i = q i /k B T with q i being the charge on species i particles and q K+1 = 0, k B is the Boltzmann constant, T is an absolute temperature, is an average volume. The following inequalities imply that the distributions are of Fermi-like type [87] C i (r) < lim i.e., C i (r) cannot exceed the maximum value 1/v 2 i or 1/v i for any arbitrary (or even infinite) potential φ(r) at any location r in the domain Ω s , where i = 1, · · ·, K + 1, The classical Boltzmann distribution appears if all particles are treated as volumeless points, i.e., v i = 0 and Γ(r) = Γ B = 1. The classical Boltzmann distribution may produce an infinite concentration C i (r) → ∞ in crowded conditions when −β i φ(r) → ∞, close to charged surfaces for example, which is physically impossible [64][65][66]. This is a major, even crippling deficiency of PB theory for modeling a system with strong local electric fields or interactions. The difficulty in the application of classical Boltzmann distributions to saturating systems has been avoided in the physiological literature (apparently starting with Hodgkin, Huxley, and Katz [88]) by redefining the Boltzmann distribution to deal with systems that can only exist in two states. This redefinition has been vital to physiological research and is used in hundreds of papers [89,90], but confusion results when the physiologists' saturating two-state Boltzmann is not kept distinct from the unsaturating Boltzmann distribution of statistical mechanics [91].
It should be clearly understood that as beautiful as is Hodgkin's derivation it begs the question of what physics creates and maintains two states. Indeed, it is not clear how one can define the word state in a usefully unique way in a protein of enormous molecular weight with motions covering the scale from femtoseconds to seconds.
The steric potential S trc (r) in Equation (4) first introduced in [64] is an entropic measure of crowding or emptiness of particles at r. If φ(r) = 0 and C i (r) = C B i then S trc (r) = 0. The factor v i /v 0 shows that the steric energy −v i v 0 S trc (r)k B T of a type i particle at r depends not only on the steric potential S trc (r) but also on its volume v i similar to the electric energy β i φ(r)k B T depending on both φ(r) and q i [72]. The steric potential S trc (r) is especially relevant to determining selectivity of specific ions by certain biological ion channels [65,66,68,70,72].
In this mean-field Fermi distribution, it is impossible for a volume v i to be completely filled with particles, i.e., it is impossible to have v i C i (r) = 1 (and thus Γ(r) = 0) since that would make S trc (r) = −∞ and hence C i (r) = 0, a contradiction. Therefore, we must include the voids as a separate species if we treat ions and water having volumes in a model for which C i (r) < 1/v i and Γ(r) = 0 for all i = 1, · · ·, K + 1 and r ∈ Ω s . This is a critical property distinguishing our theory from others that do not consider water as a molecule with volume and so do not have to consider voids. We shall elaborate this property in Section 4.
Our theory is consistent with the classical theory of van der Waals in molecular physics, which describes nonbond interactions between any pair of atoms as a distance-dependent potential such as the Lennard-Jones (L-J) potential that cannot have zero distance between the pair [85,86]. Indeed, the steric potential S trc (r) can be written as a function of the volume of all molecular species (of course, including water as well as ions). Classical extensions of van der Waals theories often use this variable, but seem not to mention the existence or importance of voids.
The steric potential S trc (r) lumps all van der Waals potential energies of paired particles in a mean-field sense. It is an approximation of L-J potentials that describe local variations of L-J distances (and thus empty voids) between any pair of particles. L-J potentials are highly oscillatory and extremely expensive and unstable to compute numerically [66]. Calculations that involve L-J potentials [92][93][94][95][96][97][98], or even truncated versions of L-J potentials [99][100][101] must be extensively checked to be sure that results do not depend on irrelevant parameters. Any description that uses L-J potentials has a serious problem specifying the combining rule. The details of the combining rule directly change predictions of effects of different ions (selectivity) and so predictions depend on the reliability of data that determines the combining rule and its parameters.
The steric potential does not require combining rules. Since we consider specific sizes of ions and water with voids, the steric potential is valid on the atomic scale of L-J potentials. It is also consistent with that on the macroscopic scale of continuum models as shown in Sections 6 and 8.
To our surprise during the writing of this article, we found Equation (2) in Bikerman's 1942 paper [42] is exactly the same as Equation (4) for a special case of binary ionic liquids with the identical steric energies of different-sized ions, i.e., the factor v i /v 0 = 1 in (4). The steric potential in Equation (4) is however not explicitly expressed in Bikerman's paper. Therefore, Bikerman's concentration function is a Fermi distribution, a generic term used in statistical mechanics. We do NOT use exactly the Fermi distribution as Fermi derived in 1926 for identical particles now called fermions in quantum mechanics. So it is both more precise and historically correct to use the name "Poisson-Bikerman" equation for finite-sized ions as a generalization of the Poisson-Boltzmann equation for volumeless ions in electrochemical and bioelectric systems.
As noted by Bazant et al. in their review paper [14], Bikerman's paper has been poorly cited in the literature until recently. In our intensive and extensive study of the literature since 2013 [64], we have never found any paper specifically using Bikerman's formula as Equation (4), although of course there may be an instance we have not found. We thus now change the term "Poisson-Fermi" used in our earlier papers to "Poisson-Bikerman" in honor of Bikerman's brilliant work. We present here mathematical as well as physical justifications of a very general treatment of different-sized ions and water molecules in the mean-field framework based on Bikerman's pioneer work.

Fourth-Order Poisson-Bikerman Equation and Correlations
Electrolytes have been treated mostly in the tradition of physical chemistry of isolated systems that proved so remarkably successful in understanding the properties of ideal gases in atomic detail, long before the theory of partial differential equations, let alone numerical computing was developed. Most applications of ionic solutions however involve systems that are not at all isolated. Rather, most practical systems include electrodes to deliver current and control potential, and reservoirs to manipulate the concentrations and types of ions in the solution. Indeed, all biology occurs in ionic solutions and nearly all of biology involves large flows. It is necessary then to extend classical approaches so they deal with external electric fields and other boundary conditions and allow flow so the theory can give useful results that are applicable to most actual systems. When the electrolyte system in Ω s is subject to external fields such as applied voltages, surface charges, and concentration gradients on the boundary ∂Ω s , the electric field E(r) of the system, the displacement field D(r) of free ions, and the polarization field P(r) of water are generated at all r in Ω s . In Maxwell's theory [102,103], these fields form a constitutive relation D(r) = 0 E(r) + P(r) (9) and the displacement field satisfies where 0 is the vacuum permittivity, ρ ion (r) is the charge density of ions, and C i (r) are the concentrations defined in (4). See [104] for a modern formulation of Maxwell's theory applicable wherever the Bohm version of quantum mechanics applies [105,106]. The electric field E(r) is thus screened by water (Bjerrum screening) and ions (Debye screening) in a correlated manner that is usually characterized by a correlation length l c [77,78,107]. The screened force between two charges in ionic solutions (at r and r in Ω s ) has been studied extensively in classical field theory and is often described by the van der Waals potential kernel [71,72,84,107,108] that satisfies the Laplace-Poisson equation [108] − ∆W(r − r ) in the whole space R 3 , where ∆ = ∇ · ∇ = ∇ 2 is the Laplace operator with respect to r and δ(r − r ) is the Dirac delta function at r . The potential φ(r) defined in describes an electric potential of free ions [72,107] that are correlated only by the mean electric field according to the Poisson equation a second-order partial differential equation, where s = w 0 and w is the dielectric constant of water. This potential does not account for correlation energies between individual ions or between ion and polarized water in high field or crowded conditions under which the size and valence of ions and the polarization of water play significant roles [17,[65][66][67][68]77,78,107]. The correlations implicit in Maxwell's equations are of the mean-field and can be summarized by the statement that current is conserved perfectly and universally on all scales that the Maxwell equations are valid, where current includes the term 0 ∂E(r,t) ∂t . This term allows the Maxwell equations to describe the propagation of light through a vacuum, and it allows charge to be relativistically invariant, i.e., independent of velocity unlike mass, length, and time all of which vary dramatically as velocities approach the speed of light [104,106].
We introduce the correlated electric potential (15) in [72] as a convolution of the displacement potential φ(r ) with W(r − r ) to deal with the correlation and polarization effects in electrolyte solutions. However, it would be too expensive to calculate φ(r) using (15). Multiplying (12) by φ(r ) and then integrating over R 3 with respect to r [71], we obtain a Laplace-Poisson equation [107,108] that satisfies (15) in the whole unbounded space R 3 with the boundary conditions φ(r) = φ(r) = 0 at infinity. From (14) and (16), we obtain the 4th-order Poisson-Bikerman equation a PDE that is an approximation of (16) in a bounded domain Ω s ⊂ R 3 with suitable boundary conditions (see below) of φ(r) on ∂Ω s . We can thus use (9) to find the polarization field with E(r) = −∇φ(r). If l c = 0, we recover the standard Poisson Equation (14) and the standard polarization P = 0 ( w − 1)E with the electric susceptibility w − 1 (and thus the dielectric constant w ) if water is treated as a time independent, isotropic, and linear dielectric medium [103]. In this case, the field relation D = w 0 E with the scalar constant permittivity s 0 is an approximation of the exact relation (9) due to the simplification of the dielectric responses of the medium material to the electric field E [109][110][111].
The exponential van der Waals potential W(r − r ) = e −|r−r |/lc |r−r |/l c [84] is called the Yukawa [112] potential in [71,72] and usually in physics, which is an anachronism [108,113]. Van der Waals derived this potential in his theory of capillarity based on the proposition that the intermolecular potential of liquids and gases is shorter-ranged, but much stronger than Coulomb's electric potential [108]. Ornstein and Zernike (OZ) introduced short-(direct) and long-ranged (indirect) correlation functions in their critical point theory [114]. There are three important properties of the van der Waals potential: (i) it satisfies the Laplace-Poisson Equation (12), (ii) it generates the same functional form for shortand long-ranged correlations in the OZ theory, and (iii) it solves van der Waals's problem for the intermolecular potential [108].
Therefore, the potential φ(r) in (15) includes correlation energies of ion-ion and ion-water interactions in short as well as long ranges in our system. The correlation length l c can be derived from the OZ equation, see Equation (13) in [108], but the derivation is not very useful. The correlation length becomes an unknown functional of ρ ion (r) in (10) and the OZ direct correlation function, and is hence usually chosen as an empirical parameter to fit experimental, molecular dynamics (MD), or Monte Carlo (MC) data [14,[64][65][66][67][68][69][70][71][72][73]75,77,78,107]. It seems clear that it would be useful to have a theory that showed the dependence of correlation length on ion composition and concentration, and other parameters.
There are several approaches to fourth-order Poisson-Boltzmann equations for modeling ion-ion and ion-water correlations from different perspectives of physics [71,77,78,115,116]. In [77], a decomposed kernel acts on a charge density of counterions in a binary liquid without volumes and water (ion-ion correlations) in contrast to the potential φ(r) in (15) that is generated by different-sized ions and water with voids in (14) (ion-ion and ion-water correlations in a multi-component aqueous electrolyte). The kernel consists of short-range (of van der Waals type) and long-range components from a decomposition of Coulomb's interactions. In [78], the kernel is a general nonlocal kernel that acts on a charge density of equal-sized ions in a binary liquid without water (ion-ion correlations). The kernel is a series expansion of the gradient operator ∇ and thus can yield not only a fourth-order PB but even higher-order PDEs. The fourth-order PB is the first-order approximation of the energy expansion that converges only with small wavenumbers k in Fourier frequency domain for the dielectric response of ionic liquids [78].
Derived from the framework of nonlocal electrostatics for modeling the dielectric properties of water in [107], the kernel acting on φ(r) in [71] (ion-ion and ion-water correlations) consists of a van der Waals function and the Dirac delta function that correspond to the limiting cases k = 0 and k = ∞, respectively. In [115], a system of three PDEs derived from electrostatics and thermodynamic pressure has electric potential and concentration gradients of equal-sized cations and anions in a binary fluid as three unknown functions. Linearization and simplification of the nonlinear system can yield a linear fourth-order PB (ion-ion correlations). In [116], the fourth-order PB is derived from a free energy functional that models ion-ion correlations in a binary liquid using volume-fraction functions of equal-sized cations and anions with two additional parameters associated with the interaction energies of these two functions and their gradients.
The dielectric operator s l 2 c ∆ − 1 in (17) describes changes in dielectric response of water with salt concentrations (ion-water correlations), ion-ion correlations, and water polarizations all via the mean-field charge density function ρ ion (r) provided that we can solve (4) and (17) for a consistent potential function φ(r). Therefore, the operator (a mapping) depends not only on ion and water concentrations (C B i for all arbitrary species i = 1, · · ·, K + 1 of particles with any arbitrary shapes and volumes) but also on the location r and the voids at r. The operator thus produces a dielectric function (r,C B i ) as an output from the solution φ(r) that satisfies the 4PBik (17) that saturates as a function of concentration (4), as we shall repeatedly emphasize. This dielectric function (r,C B i ) is not an additional model for (r), (k), or (C B i ) as it often is in other models in the literature [62,63,[117][118][119][120][121][122][123][124][125][126][127]. Here the dielectric function is an output, as we have stated.
The 4PBik Equation (17) with (4) is a very general model using only one extra parameter l c in the fourth-order operator to include many physical properties ignored by the classical Poisson-Boltzmann equation. We shall illustrate these properties of our model in Section 8.

Generalized Gibbs Free Energy Functional
To generalize the Gibbs free energy functional for Boltzmann distributions that satisfy the classical Poisson-Boltzmann equation [3,128,129], we introduce a functional in [72] for saturating Fermi distributions (4) that satisfy the 4 th -order Poisson-Bikerman Equation (17) where F el (C) is an electrostatic functional, F en (C) is an entropy functional, C = (C 1 (r), C 2 (r), · · ·, C K+1 (r)), and L −1 is the inverse of the self-adjoint positive linear operator (17), i.e., Lφ(r) = ρ ion (r). C is a 'concentration vector' that specifies the number density, i.e., concentration of each species in the ionic solution, including water. C plays a central role in any theory of ionic solutions because it specifies the main property of a solution, namely its composition.
Taking the variations of F(C) at C i (r), we have that yields the saturating Fermi distributions in (4) for all i = 1, · · ·, K + 1. Moreover, we have implying that the saturating Fermi distribution vector C is a unique minimizer of the functional F(C).
The Gibbs-Bikerman free energy functional F(C) has two important properties. First, its electrostatic part F el (C) is defined in terms of the composition vector C only. It depends only on concentrations and nothing else. If an electrostatic functional F el ( φ(r)) is defined in terms of ∇ φ(r) 2 for the PB equation [64,78,124,[130][131][132][133][134][135][136][137], the corresponding concentration vector C and the potential φ(r) do not minimize the corresponding functional F( C, φ(r)) [128,129], i.e., F is not a Gibbs free energy functional [3,128]. Second, the limit of its entropic part exists (F en converges) when the volume v i tends to zero for all i = 1, · · ·, K + 1. This implies that all ionic species have Boltzmann distributions C 0 is a constant, and the void fraction Γ(r) = Γ B = 1 since all particles are volumeless in PB theory. Therefore, the 4PBik model (4) and (17) is physically and mathematically consistent with the classical PB model in the limiting case when we ignore the steric (v i = 0) and correlation (l c = 0) effects.
There are many shortcomings of the lattice approach [138] frequently used to account for steric effects in lattice-based PB models [14,61,64,78,124,129,[133][134][135][136][139][140][141]. For example, (i) it assumes equal-sized ions and thus cannot distinguish non-uniform particles as in (1), (ii) its effective ion size needs to be unrealistically large to fit data [14], (iii) its correction over Boltzmann's point charge approach appears only at high surface charges [125], (iv) its pressure term diverges very weakly (is greatly underestimated) at close packing [142], and (v) its entropy functional may diverge as the volume of ions tends to zero, i.e., the corresponding lattice-based PB model is not physically and mathematically consistent with the classical PB model in the limiting case [66].
The importance of the restriction in Point (i) is hard to overstate. Almost all the interesting properties of ionic solutions arise because of their selectivity (as it is called in biology) or specificity between species, and those different properties arise in large measure because of the different diameters of the ions. The equal diameter case is dull and degenerate. Point (v) is a critical problem that is closely related to Points (ii)-(iv). The divergence is obvious for an entropy term F en in Equation (2) in [133] as which also appears in [61,64,78,124,129,[133][134][135][136][139][140][141]. It is impossible to derive Boltzmann [129], for example. In fact, the assumption (2.6), i.e., v C i (r) > 0, actually forbids us from taking v to the limit zero.
Our derivation of F en (C) does not employ any lattice models but simply uses the exact volume Equation (1). Our theory should not be classified then as a lattice model as sometimes is the case, at least in informal discussions. The void function Γ(r) is an analytical generalization of the void fraction 1 − Φ in (20) in [14] with all volume parameters v i (including the bulk fraction Γ B ) being physical instead of empirical as Φ. The excess chemical potential in [14] is −k B T ln(1 − Φ) whereas ours is F en (C) in (21).
These expressions are different in important respects. Our model is not a lattice-based model because its differences are crucial both mathematically and physically. Indeed, the lattice-based model is in a certain sense internally inconsistent with classical statistical mechanics since a fundamental result of classical statistical mechanics v C i (r) > 0 prevents the model from satisfying the classical imperative of the Boltzmann distribution in the limit of zero v.
The Langmuir-type distribution of different-sized ions (without water) proposed in [125] also reduces to a Boltzmann distribution as v j → 0, ∀j, where C max j = p/v j and p ≤ 1 is a packing parameter. This distribution saturates and thus is of Fermi type, i.e., (4). Our distribution in (4) does not need any packing parameters and satisfies v i C i (r) < 1.

Poisson-Nernst-Planck-Bikerman Model of Saturating Phenomena
For nonequilibrium systems, we can also generalize the classical Poisson-Nernst-Planck model [38,39,43,143,144] to the Poisson-Nernst-Planck-Bikerman model by coupling the flux density equation of each particle species i = 1, · · ·, K + 1 (including water) to the 4PBik Equation (17), where the flux density is defined as D i is the diffusion coefficient, and the time variable t is added to describe the dynamics of electric φ(r, t) and steric S trc (r, t) potentials. The flux Equation (27) is called the Nernst-Planck-Bikerman equation because the steric potential S trc (r, t) is introduced into the classical NP equation so it can deal with saturating phenomena including those that arise from the unequal volumes of ions and the finite volume of molecular water. The PNPB model can be extended to include hydrodynamic kinetic and potential energies in the variational treatment of energy processes (i.e., EnVarA) by Hamilton's least action and Rayleigh's dissipation principles [145,146]. We shall however consider this as a topic for future work.
At equilibrium, the net flow of each particle species is a zero vector, i.e., J i (r) = 0 (in a steady state), which implies that where the constant c i = C B i for φ(r) = S trc (r) = 0. Therefore, (29) = (4), i.e., the NPB Equation (27) reduces to the saturating Fermi distribution (4) as the classical NP equation reduces to the Boltzmann distribution at equilibrium.
The gradient of the steric potential ∇S trc (r, t) in (28) represents an entropic force of vacancies exerted on particles. The negative sign in −C i (r, t)∇S trc (r, t) means that the steric force ∇S trc (r, t) is in the opposite direction to the diffusion force ∇C i (r, t).
Larger S trc (r, t) = ln Γ(r,t) Γ B implies lower pressure because the ions occupy more space (less crowded) as implied by the numerator Γ(r, t). The larger the S trc (r, t) the lower pressure at the location r, the more the entropic force (the higher pressure) pushes particles to r from neighboring locations. The steric force is the opposite of the diffusion force ∇C i (r, t) that pushes particles away from r if the concentration at r is larger than that at neighboring locations.
Moreover, the Nernst-Einstein relationship between diffusion and mobility [9] implies that the steric flux where the mobility coefficient µ i of an ion depends on its size v i in addition to its charge q i .
Therefore, the gradients of electric and steric potentials (∇φ(r, t) and ∇S trc (r, t)) describe the charge/space competition mechanism of particles in a crowded region within a mean-field framework. Since S trc (r, t) describes the dynamics of void movements, the dynamic crowdedness (pressure) of the flow system can also be quantified. A large amount of experimental data exists concerning the dependence of diffusion coefficient on the concentration and size of solutes. Comparing our model with this data is an important topic of future work.
The motion of water molecules, i.e., the osmosis of water [147,148] is directly controlled by the steric potential in our model and their distributions are expressed by Nevertheless, this motion is still implicitly changed by the electric potential φ(r, t) via the correlated motion of ions described by other C j (r, t) in the void fraction function Γ(r, t) and hence in the charge density ρ ion (r, t) in (17).
In summary, the PNPB model accounts for (i) the steric (pressure) effect of ions and water molecules, (ii) the correlation effect of crowded ions, (iii) the screening (polarization) effect of polar water, and (iv) the charge/space competition effect of ions and water molecules of different sizes and valences. These effects are all closely related to the interstitial voids between particles and described by two additional terms, namely the correlation length and the steric potential. The steric potential is most naturally written as a function of the volume of voids, but it can also be written as a function of the total volume of all molecules, including water and ions.

Generalized Debye-Hückel Theory
Thermodynamic modeling is of fundamental importance in the study of chemical and biological systems [1,6,9,[11][12][13]16,32]. Since Debye and Hückel (DH) proposed their theory in 1923 [149] and Hückel extended it to include Born energy effects in 1925 [150], a great variety of extended DH models (equations of state) have been developed for modeling aqueous or mixed-solvent solutions over wide ranges of composition, temperature, and pressure [6,19,[151][152][153][154][155]. Despite these intense efforts, robust thermodynamic modeling of electrolyte solutions still presents a difficult challenge for extended DH models due to an enormous number of parameters that need to be adjusted carefully and often subjectively [19,[152][153][154]156].
It is indeed a frustrating despair (the word frustration on p. 11 in [16] and the word despair on p. 301 in [1]) that about 22,000 parameters [19] need to be extracted from the available experimental data for one temperature for combinatorial solutions of the most important 28 cations and 16 anions in salt chemistry by the Pitzer model [6], which is the most widely used DH model with unmatched precision for modeling electrolyte solutions [153]. The JESS (joint expert speciation system) is the world's largest system of thermodynamic information relating to electrolytes, reactions in aqueous media, and hydrocarbon phase equilibria [157]. The total number of Pitzer's fitting parameters in JESS is 95 [158].
We briefly outline the derivation of a generalized DH equation of state and refer to [76] for more details. The activity coefficient γ i of an ion of species i in an aqueous electrolyte solution with a total of K species of ions describes deviation of the chemical potential of the ion from ideality (γ i = 1) [11]. The excess chemical potential µ ex i = k B T ln γ i can be calculated by [69,167] where q i is the charge of the hydrated ion (also denoted by i), φ(r) is a reaction potential [167] function of spatial variable r in the domain Ω = Ω i ∪ Ω sh ∪ Ω s shown in Figure 1, Ω i is the spherical domain occupied by the ion i, Ω sh is the hydration shell domain of the ion, Ω s is the rest of solvent domain, 0 denotes the center (set to the origin) of the ion, and φ 0 (r) is a potential function when the solvent domain Ω s does not contain any ions at all with pure water only, i.e., when the solution is ideal. The radii of Ω i and the outer boundary of Ω sh are denoted by R Born i (ionic cavity radius [160]) and R sh i , respectively. The potential function φ(r) can be found by solving the 4PBik Equation (17) and the Laplace equation [69,73] ∆φ where s is defined in Ω sh ∪ Ω s , the correlation length l c = √ l B l D /48 is a density-density correlation length independent of specific ionic radius [168], l B and l D are the Bjerrum and Debye lengths, respectively, the concentration C k (r) function (4) is defined in Ω for all k = 1, · · ·, K + 1 in molarity (M), and v k = 4πa 3 k /3 with radius a k . Since the steric potential takes particle volumes and voids into account, the shell volume V sh of the shell domain Ω sh can be determined by the steric potential 69,73], where the occupant (coordination) number O w i of water molecules is given by experimental data [166]. The shell radius R sh i is thus determined and depends not only on O w i but also on the bulk void fraction Γ B , namely on all salt and water bulk concentrations (C B k ) [69,73]. We reduce the complexity of higher-order approximations, and make them easier to implement by transforming the fourth-order PDE (17) to the following two second-order PDEs [64] where the extra unknown function ψ(r) is a density-like function as seen from (33) where ∂ denotes the boundary of a domain, the jump function [φ(r)] = lim r sh →r φ(r sh ) − lim r i →r φ(r i ) at r ∈ ∂Ω i with r sh ∈ Ω sh and r i ∈ Ω i , (r) = s in Ω sh and (r) = ion 0 in Ω i , ion is a dielectric constant in Ω i , n is an outward normal unit vector at r ∈ ∂Ω i , and φ * (r) = q i /(4π i |r − 0|). Equation (32) avoids large errors in a direct approximation of the delta function δ(r − 0) in the singular charge q i δ(r − 0) of the solvated ion at the origin 0 by transforming the singular charge to the Green's function φ * (r) on ∂Ω i in (39) as an approximation source of the electric field produced by the solvated ion [169,170]. For simplicity, we consider a general binary (K = 2) electrolyte C z 2 A z 1 with the valences of the cation C z 1 + and anion A z 2 − being z 1 and z 2 , respectively. The first-order Taylor approximation of the charge density functional ρ ion (φ(r)) in (17) with respect to the electric potential φ(r) yields where which is a quantity corresponding to a linearization of the steric potential S trc (r) [76]. Consequently, we obtain a generalized Debye length that reduces to the original Debye length l D [11] if v 1 = v 2 = 0 (two ionic species with equal radius and thus Λ = 0) or v 1 = v 2 = v 3 = 0 (all particles treated as volumeless points in standard texts for PB [11]). The nonlinear value of Λ = 0 for v 1 = v 2 = 0 can be obtained by Newton's method [76]. Equation (33) is a second-order PDE that requires two boundary conditions like (35) and (36) for a unique solution ψ(r). Since ψ(r) = s ∇ 2 φ(r) = −ρ(r) ≈ s κ 2 φ(r) if l c = 0, Equation (36) is a simplified (approximate) boundary condition for ψ(r) on ∂Ω sh ∩ ∂Ω s without involving higher-order derivatives of ψ(r) (or the third-order derivative of φ(r)). The approximations in (36) and (40) do not significantly affect our generalized DH model's ability to fit activity data. However, these assumptions should be carefully scrutinized in other applications such as highly charged surfaces. Bazant et al. have recently developed more consistent and general boundary conditions for their fourth-order model by enforcing continuity of the Maxwell stress at a charged interface [171,172].
In [76], we analytically solve the linear 4PBik PDEs (32), (33), and (34) with (40) in a similar way as Debye and Hückel solved the linear PB equation for a spherically symmetric system. However, the spherical domain shown in Figure 1 and the boundary and interface conditions in (35)- (39) are different from those of the standard method for the linear PB equation in physical chemistry texts [11]. The analysis consists of the following steps: (i) The nonlinear term ρ ion (r) in (33) The analytical potential function that we found [76] is in Ω s , where [76]. The linearized 4PBik potential φ 4PBik (r) reduces to the linearized PB potential φ PB (r) = q i e −r/l D /(4π s r) as in standard texts (e.g., Equation (7.46) in [11]) by taking lim l c →0 φ 4PBik (r) with v k = 0 for all k, R sh i = 0, and r > 0 [76]. As discussed in [173], since the solvation free energy of an ion i varies with salt concentrations, the Born energy q 2 where C B i = C B i /M is a dimensionless bulk concentration and α i 1 , α i 2 , and α i 3 are parameters for modifying the experimental Born radius R 0 i to fit experimental activity coefficient γ i that changes with the bulk concentration C B i of the ion. The Born radii R 0 i given below are from [173] obtained from the experimental hydration Helmholtz free energies of those ions given in [12]. The three parameters in (44) have physical or mathematical meanings unlike many parameters in the Pitzer model [19,153,156]. The first parameter α i 1 adjusts R 0 i and accounts for the real thickness of the ionic atmosphere (Debye length), which is proportional to the square root of the ionic strength in the DH theory [11]. The second α i 2 and third α i 3 parameters are adjustments in the next orders of approximation beyond the DH treatment of ionic atmosphere [73].
The potential value φ 0 (0) (31) and (42), we thus have a generalized activity coefficient γ 4PBik if R Born i = R 0 i (without considering Born energy effects), R sh i = R i (an effective ionic radius [149]), l D4PBik = l D (no steric effect), and l c = 0 (no correlation effect). The reduction shown in [76] is by taking the limit of the last term in (45) as l c → 0, i.e., lim l c →0 Hückel soon realized that the DH formula (46) failed to fit experimental data at high ionic strengths and modified it in 1925 [150] by adding one more parameter η 1 to become (see Equation (7.115) in [11]) where η 0 (an approximation of R sh i ) and η 1 account for the distance of closest approach to the ion i and the salting-out effect (an approximation of the Born energy), respectively [11], where I = 1 2 ∑ i C B i z 2 i is the ionic strength of the solution. Consequently, a variety of extended DH models γ DHBx i [153,174] in the form similar to have been proposed in the literature to express other thermodynamic properties such as temperature and pressure by a power expansion of I with more and more parameters η k that can increase combinatorially with various composition, temperature, and pressure to a frustrating amount [1,16,19]. Please note that η k may also depend on ionic strength I in a complicated way, see e.g., Equation (2) in [153]. Many expressions of those parameters are rather long and tedious and do not have clear physical meaning [19,153,156]. The Davies equation [175] is a special form of (47) with a linear term in I. The R Born i term in (45) differs significantly from the last term in (48) as they are the inverse of each other in terms of I and parameters, i.e., I, α 1 , α 2 , and α 3 are in the denominator in (45) whereas I and η k are in the numerator in (48). This implies that γ 4PBik i and γ DHBx i vary oppositely with I. Consequently, the values of α 1 , α 2 , and α 3 are totally different from those of η k when we use γ 4PBik i and γ DHBx i to fit experimental activity coefficients with I varying from low to high values [76]. This may explain why the empirical nature of extended DH models requires a great deal of effort to extract parameters (without physical hints) from existent thermodynamic databases by regression analysis [19,153].

Numerical Methods
Numerical simulations are indispensable to study chemical, physical, and mathematical properties of biological and chemical systems in realistic applications, especially with experimental details at atomic scale such as ion channels in the Protein Data Bank (PDB) [57]. Continuum PDE models have substantial advantages over Monte Carlo, Brownian dynamics (BD), or molecular dynamics in physical insights and computational efficiency that are of great importance in studying a range of conditions and concentrations especially for large nonequilibrium or inhomogeneous systems, as are present in experiments and in life itself [10,17,21,95,121,[176][177][178][179][180][181][182][183][184][185].
The literature on numerical methods for solving PB and PNP models is vast [64,68,74]. We summarize here some important features of the methods proposed in [64,68,74] for Poisson-Bikerman and Poisson-Nernst-Planck-Bikerman models, which may be useful for workers in numerical analysis and coding practice. Since PNPB including 4PBik is highly nonlinear and the geometry of protein structures is very complex, we emphasize two different types of methods, namely nonlinear iterative methods and discretization methods for these two problems as follows.

Nonlinear Iterative Methods
For the PNPB system of K + 1 NP equations in (27), Laplace Equation (32), and two 4PBik equations in (33) and (34), the total number of second-order PDEs that we need to solve is K + 4. These PDEs are coupled together and highly nonlinear except (32). Numerically solving this kind of nonlinear systems is not straightforward [64,68,74]. We use the following algorithm to explain essential procedures for solving the steady-state PNPB system, where Ω m denotes the biomolecular domain that contains a total of Q fixed atomic charges q j located at r j in a channel protein as shown in Figure 2L for the gramicidin A channel downloaded from PDB with Q = 554, for example, ∂Ω m denotes the molecular surface of the protein and the membrane lipids through which the protein crosses as shown in Figure 2R, and Ω s is the solvent domain consisting of the channel pore and the extracellular and intracellular baths for mobile ions and water molecules.
in Ω s with φ New = V on ∂Ω and the same jump condition in Step 2, where ρ (φ) is the derivative of ρ(φ) with respect to φ. 6.

Solve 4PBik1 Equation for Ψ New as in
Step 4 with C New i in place of C Old i . 9.
Solve 4PBik2 Equation for φ New as in Step 5.
Linearizing the nonlinear 4PBik (17) yields two second-order linear 4PBik1 and 4PBik2 in Steps 4 and 5 that differ from the nonlinear (33) and (34). Newton's iterative Steps 4-6 for solving 4PBik1 and 4PBik2 dictates convergence that also depends on various mappings from an old solution φ Old to a new solution φ New . This algorithm uses two relaxation and three continuation mappings for which we need to carefully tune two relaxation parameters ω 4PBik and ω PNPB and three continuation parameters λ c (related to correlation effects), λ s (steric effects), and ∆V (incremental voltage for applied voltage). For example, the parameter λ s in Γ Old (r) = 1 − ∑ K+1 j=1 λ s v j C Old j (r) can be chosen as λ s = k∆λ, k = 0, 1, 2, · · ·, 1 ∆λ , an incremental continuation from 0 (no steric effects) to 1 (fully steric effects) with a tuning stepping length ∆λ. The algorithm can fail to converge if we choose ∆λ = 1 (without continuation) for some simulation cases, since we may have Γ Old (r) < 0 resulting in numerically at some r where the potential φ Old (r) is large.

Discretization Methods
All PDEs in Steps 1, 2, 4, 5, 8, and 9 are of Poisson type −∇ 2 φ(r) = f (r). We use the central finite difference (FD) method [64] to discretize it at all grid points r ijk = (x i , y j , z k ) in a domain, where φ ijk ≈ φ(x i , y j , z k ), f ijk = f (x i , y j , z k ), and ∆x, ∆y, and ∆z are mesh sizes on the three axes from a uniform partition ∆x = ∆y = ∆z = h. The domains in Steps 1 and 2 are Ω m and Ω s , respectively. The discretization leads to a sparse matrix system A − → φ = − → f with the compressed bandwidth of the matrix A being 7, where the matrix size can be millions for sufficiently small h to obtain sufficiently accurate φ ijk .
The matrix system consists of four subsystems, two by the FD method (49) in Ω m and Ω s , one by another method (see below) to discretize the jump condition in Step 2 on the interface ∂Ω m between Ω s and Ω m , and one by imposing boundary conditions on ∂Ω. We need to solve the matrix system in the entire domain Ω = Ω m ∪ Ω s .
The convergence order of (49) is O(h 2 ) (optimal) in maximum error norm for sufficiently smooth function φ(r). However, this optimal order can be easily degraded to O(h 0.37 ) [186], for example, by geometric singularities if the jump condition is not properly treated. In [64], we propose the interface method where at every jump position γ ∈ ∂Ω m that is at the middle of its two neighboring grid points, i.e., and x i−1 and x i belong to different domains Ω s and Ω m . The corresponding cases in yand z-axis follow obviously in a similar way. This method yields optimal convergence [64]. Since the matrix system is usually very large in 3D simulations and we need to repeatedly solve such systems updated by nonlinear iterations as shown in the above algorithm, linear iterative methods such as the bi-conjugate gradient stabilized (bi-CG) method are used to solve the matrix system [74]. We propose two parallel algorithms (one for bi-CG and the other for nonlinear iterations) in [74] and show that parallel algorithms on GPU (graphic processing unit) over sequential algorithms on CPU (central processing unit) can achieve 22.8× and 16.9× speedups for the linear solver time and total runtime, respectively.
Discretization of Nernst-Planck equations in Step 7 is different from (49) because the standard FD method and thereby a negative (unphysical) concentration C i < 0 at x i if where Therefore, it is crucial to check whether the generalized Scharfetter-Gummel (SG) condition [68] − β∆φ i + ∆S trc i ≤ 2 (55) is satisfied by any discretization method in implementation. This condition generalizes the well-known SG stability condition in semiconductor device simulations [187,188] to include the steric potential function S trc (r). We extend the classical SG method [187] of the flux J(x) in [68] to where t i = β∆φ i − ∆S trc i and B(t) = t e t −1 is the Bernoulli function [188]. Equation (56), an exponential fitting scheme, satisfies (55) and is derived from assuming that the flux J, the local electric field −dφ dx , and the local steric field dS trc dx are all constant in the sufficiently small subinterval (x i , x i+1 ), i.e., where k = β dφ dx − dS trc dx . Solving this ordinary differential equation (ODE) with a boundary condition C i or C i+1 yields the well-known Goldman-Hodgkin-Katz flux equation in ion channels [9], which is exactly the same as that in (56) but with the subinterval (x i , x i+1 ) being replaced by the height of the entire box in Figure 2R.

The generalized Scharfetter-Gummel method for Nernst-Planck equations is thus
The SG method is optimal in the sense that it integrates the ODE (57) exactly at every grid point with a suitable boundary condition [189]. Therefore, the SG method can resolve sharp layers very accurately [189] and hence needs few grid points to obtain tolerable approximations when compared with the primitive FD method. Moreover, the exact solution of (57) for the concentration function C(x) yields an exact flux J(x). Consequently, the SG method is current preserving, which is particularly important in nonequilibrium systems, where the current is possibly the most relevant physical property of interest [190].
It is difficult to overstate the importance of the current preserving feature and it must be emphasized for workers coming from fluid mechanics that preserving current has a significance quite beyond the preserving of flux in uncharged systems. Indeed, conservation of current (defined as Maxwell did to include the displacement current of the vacuum 0 ∂E(r,t) ∂t ) is an unavoidable consequence, nearly a restatement of the Maxwell equations themselves [104,106]. The electric field is so strong that the tiniest error in preserving current, i.e., the tiniest deviation from Maxwell's equations, produces huge effects. The third paragraph of Feynman's lectures on electrodynamics makes this point unforgettable [191]. Thus, the consequences of a seemingly small error in preserving the flow of charge are dramatically larger than the consequences of the same error in preserving the flux of mass.
We have developed a C++ code for solving 4PBik and PNPB models on laptop and highperformance (with GPU) computers. For solving a 4PBik problem with a matrix system of size 4,096,000, for example, the code requires about 300 MB memory to store the compressed matrix system with double precision. It took about 2 min and 47 s on a laptop computer equipped with 1.3 GHz Intel CPU and 2 GB RAM to solve the linear system once by the successive overrelaxation method with an error tolerance of 10 −6 [64].

Applications
We have used the saturating Poisson-Nernst-Planck-Bikerman theory to study ion activities, electric double layers, and biological ion channels in the past. The theory accounts for the steric effect of ions and water molecules, the effects of ion-ion and ion-water correlations, the screening and polarization effects of polar water, and the charge/space competition effect of ions and water molecules of different sizes and valences. These effects are all closely related to the dielectric operator in (17) and the steric potential in (4) that works on both macroscopic and atomic scales. We now illustrate these properties in the following three areas using mostly experimental data to verify the theory.

Ion Activities
The curves in Figure 3 obtained by the generalized Debye-Hückel Formula (45) [75] fit well to the experimental data by Wilczek-Vera et al. [192] for single-ion activities in 8 1:1 electrolytes. There are only three fitting parameters in the formula, namely α i 1 , α i 2 and α i 3 , which we reiterate have specific physical meaning as parameters of the water shell around ions. The values of the parameters are given in Table 1 from which we observe that R Born i deviates from R 0 i slightly. For example, R Born Cl − /R 0 Cl − = 1.007∼1.044 (not shown) for Figure 3a with [LiCl] = 0∼2.5 M, since the cavity radius R Born Cl − is an atomic measure from the infinite singularity δ(r − 0) at the origin, i.e., φ 4PBik (r) and thus γ 4PBik i are very sensitive to R Born i . On the other hand, γ 4PBik i is not very sensitive to R sh i (R sh Cl − = 5.123∼5.083 Å), i.e., the fixed choice of O w i = 18 (an experimental value in [166]) for all curves is not critical but reasonable [69]. The error between the estimated O w i and its unknown true value can always be compensated by small adjustments of R Born i . The values of other symbols are a Li + = 0.6 Å, a Na + = 0.95 173], w = 78.45, ion = 1, T = 298.15 K, where a i is the (Pauling) radius of type i particle (ion) [173]. Table 1 also shows the significant order of these parameters, i.e., α i 1 > α i 2 > α i 3 in general cases. Please note that the three sets of the values of α Na + 1 , α Na + 2 , and α Na + 3 for the same Na + in three different salts NaCl, NaBr, and NaF are different because the diameters of the anions are different. Figure 4 shows single-ion activities in 6 2:1 electrolytes by experiments [192] and 4PBik, where the significant order (not shown) of three fitting parameters is similar to that in Table 1.  The electric potential and other physical properties of ionic activity can be studied in detail according to the partitioned domain in Figure 1  Li + (0) and φ 4PBik K + (0) is due to the sizes of Li + and K + not Br − as it is the same for both solutions. This example clearly shows the atomic properties of 4PBik theory in the ion Ω i and shell Ω sh domains and the continuum properties in the solvent domain Ω s . The Born radius R Born i in (42) determined by (44) changes with (i) ion-water interactions in Ω i ∪ Ω sh and (ii) ion-ion interactions in Ω i ∪ Ω s via φ 4PBik (r) in (42) that is self-consistently determined by the interface conditions in (35)- (39) and by (iii) multi-salt [73,76] concentrations in Ω s , (iv) the screening effects of water in Ω sh and ions and water in Ω s , (v) the polarization effect of water in Ω s , (vi) the correlation effect between ions in Ω s , (vii) the steric effects of all ions and water in the entire domain Ω = Ω i ∪ Ω sh ∪ Ω s , (viii) temperatures [73,76], and (ix) pressures [73,76]. The generalized Debye-Hückel formula (45) includes all these 9 physical properties with only 3 fitting parameters. However, we look forward to the day when we can derive the three fitting parameters for particular types of ions, from independently determined experimental data.

Electric Double Layers
We consider a charged surface in contact with a 0.1 M 1:4 aqueous electrolyte, where the charge density is σ = 1e/(50 Å 2 ), the radius of both cations and anions is a = 4.65 Å (in contrast to an edge length of 7.5 Å of cubical ions in [133]), and s = 80 [72]. The multivalent ions represent large polyanions adsorbed onto a charged Langmuir monolayer in experiments [133]. We solve (33) and (34) using (49) in the rectangular box Ω = Ω s = (x, y, z) : 0 ≤ x ≤ 40, − 5 ≤ y ≤ 5, − 5 ≤ z ≤ 5 Å such that φ(r) ≈ 0 within the accuracy to 10 −4 near and on the surface x = 40 Å. The boundary conditions on the surface and its adjacent four planes are − s ∇φ · n = σ with n = −1, 0, 0 and − s ∇φ · n = 0 with n defined similarly, respectively. The classical PB model (with a = a H 2 O = l c = 0, i.e., no size, void, and correlation effects) produces unphysically high concentrations of anions (A 4− ) near the surface as shown by the dashed curve in Figure 6L. The dotted curve in Figure 6L is similar to that of the modified PB in [133] and is obtained by the 4PBik model with l c = 0 (no correlations), V K+2 = 0 (no voids), and a H 2 O = 0 (water is volumeless as in [133] and hence i is the bulk water volume fraction). The voids (V K+2 = 0) and water molecules (a H 2 O = 0) have slight effects on anion concentration (because of saturation) and electric potential (because water and voids have no charges) profiles as shown by the thin solid curves in Figure 6L,R, respectively, when compared with the dotted curves. However, ion-ion correlations (with l c = 1.6a [78]) have significant effects on ion distributions as shown by the thick solid and dash-dotted curves in Figure 6L, where the saturation layer is on the order of ionic radius a and the overscreening layer [78] is about 18 Å in thickness. The saturation layer is an output (not an imposed condition) of our model unlike a Stern layer [193] imposed by most EDL models to account for size effects near charge surfaces [194][195][196]. The electric potentials φ(0) = 5.6 k B T/e at x = 0 and φ(11.5) = −1.97 k B T/e in Figure 6R obtained by 4PBik with voids and correlations deviate dramatically from those by previous models for nearly 100% at x = 0 (in the saturation layer) and 70% at x = 11.5 Å (in the screening layer) when compared with the maximum potential φ(0) = 2.82 k B T/e of previous models. The 4PBik potential depth φ(11.5) = −1.97 k B T/e of the overscreening layer is very sensitive the size a of ions and tends to zero as a → 0.

Biological Ion Channels
Biological ion channels are a particularly appropriate test of a model of concentrated ionic solutions.
The data available for tens to hundreds of different types of channels and transporters is breathtaking: it is often accurate to a few per cent (because signal to noise ratios are so large and biological variation hardly exists for channels of known amino acid sequence, which means nearly every channel presently). The data is always nonequilibrium, i.e., current voltage relations in a wide range of solutions of different composition and concentration, or (limiting zero voltage) conductance in those solutions. Indeed, many of the channels do not function if concentrations are equal on both sides and the electrical potential is zero. They are said to inactivate.
The data is often available for single channels recorded individually in patch clamp or bilayer configuration. Data is available for a range of divalent (usually calcium ion) concentrations because calcium concentration is often a controller of channel, transporter, and biological activity in the same sense that a gas pedal is the controller of the speed of a car. The structure of the ion channel or transporter is often known in breathtaking detail. The word 'breathtaking' is appropriate because similar structures are rarely if ever known of strictly physical systems. The structure and the structure of the permanent and polarization charge of the channel protein (that forms the pore through which ions move) can be modified by standard methods of site directed mutagenesis, for example that are available in 'kit' form usable by most molecular biology laboratories. Thus, models can be tested from atomic detail to single-channel function to ensemble function to cellular and physiological function, even to the ultimate biological function (like the rate of the heartbeat). Few other systems allow experimental measurement at each level of the hierarchy linking the atomic composition of genes (that encode the channel's amino acid composition), to the atomic structure of the channel, right to the function of the cell. The hierarchy here reaches from 10 −11 to 10 −5 m. When the channel controls the biological function of an organ like the heart, the hierarchy reaches to 2 · 10 −1 m, in humans for example.
The biological significance of ion channels is hard to exaggerate since they play a role in organisms analogous to the role of transistors in computers. They are the device that execute most of the physical controls of current and ion movement that are then combined in a hierarchy of structures to make biological cells, tissues, and organisms, if not populations of organisms.
From a physical point of view, ion channels provide a particularly crowded environment in which the effects of the steric potential (crowding in more traditional language) and electrical potential can combine to produce striking characteristics of selectivity and rectification. Theories that do not deal explicitly with ion channel data, i.e., that do not predict current voltage relations from known structures, seem to us to be begging central PHYSICAL questions that might falsify their approach. In fact, as a matter of history it is a fact that we learned how to construct our model of bulk solutions from our earlier work on ion channels.

Gramicidin A Channel
We use the gramicidin A (GA) channel in Figure 2L to illustrate the full Poisson-Nernst-Planck-Bikerman theory consisting of Equations (4), (27), (28), (32)- (34), and conditions (35)-(39) with-steric, correlation, polarization, dielectric, charge/space competition, and nonequilibrium effects-at steady state using the algorithm and methods in Section 7 to perform numerical simulations. The union domain Ω i ∪ Ω sh in Figure 1 is replaced by the biomolecular domain Ω m in Figure 2R. Figure 7L shows I-V curves obtained by PNPB and compared with experimental data (symbols) by Cole et al. [197] with bath K + and Cl − concentrations C B = 0.1, 0.2, 0.5, 1, 2 M and membrane potentials ∆V = 0, 50, 100, 150, 200 mV. The PNPB currents in pico ampere (pA) were obtained with θ = 1/4.7 in the pore diffusion coefficients θD i from (30) for all particle species. The reduction parameter θ has been used in all previous PNP papers and is necessary for continuum is comparable to MD, BD, or experimental data [198]. The values of other model parameters are listed in Table I in [68]. We summarize the novel results of PNPB in [68] when compared with those of earlier continuum models for ion channels: (i) The pore diffusion parameter θ = 1/4.7 agrees with the range 1/3 to 1/10 obtained by many MD simulations of various channel models [199][200][201] indicating that the steric ( Figure 7R), correlation, dehydration ( Figure 8L), and dielectric ( Figure 8R) properties have made PNPB simulations closer (realistic) to MD than previous PNP for which θ differs from MD values by an order to several orders of magnitude [200]. (ii) Figures 7R and 8L,R, which are all absent in earlier work, show that these properties correlate to each other and vary with salt concentration and protein charges in a self-consistent way by PNPB. (iii) The steric potential profiles in Figure 7R clearly illustrate the charge/space competition between ions and water under dynamic and variable conditions. For example, the global minimum value in Figure 7R at r = 13.1 on the channel axis, where the channel protein is most negatively charged, is S trc ( r) = ln Γ( r) Γ B = −0.485 yielding Γ( r)/Γ B = 0.616. Namely it is 38.4% more crowded at r than in the bath and mainly occupied by K + as shown in Figures 8L and  9L. It is important to quantify voids (Γ(r) = 1 − ∑ K+1 i=1 v i C i (r)) at highly charged locations in channel proteins and many more biological, chemical, and nano systems. The charge space competition has been a central topic in the study of ion channels since at least [202][203][204][205][206]. The literature is too large to describe in detail here. Recent reviews can help [207][208][209][210]. (iv) PNPB preserves mass conservation due to void and size effects in contrast to PNP as shown in Figure 9R, where the total number of H 2 O and K + in the channel pore is 8 [211].   [211], which is conserved by PNPB but not by PNP (without steric and correlation effects).

L-Type Calcium Channel
L-type calcium channels operate very delicately in physiological and experimental conditions. They exquisitely tune their conductance from Na + -flow, to Na + -blockage, and to Ca 2+ -flow when bath Ca 2+ varies from trace to high concentrations as shown by the single-channel currents in femto ampere in Figure 10L (circle symbols) recorded by Almers and McCleskey [212], where the range of extracellular concentrations [Ca 2+ ] o is 10 8 -fold from 10 −10 to 10 −2 M. We used the Lipkind-Fozzard molecular model [213] shown in Figure 10R to perform PNPB simulations with both atomic and continuum methods (Algorithm 2 in [68]) for this model channel, where the EEEE locus (four glutamate side chains modeled by 8 O 1/2− ions shown by red spheres) forms a high-affinity Ca 2+ binding site (center violet sphere) that is essential to Ca 2+ selectivity, blockage, and permeation. Water molecules are shown in white and red. More realistic structures would be appropriate if the work were done now, but the analysis here shows the ability of PNPB to deal with experimental data using even a quite primitive model of the structure.
PNPB results (plus symbols) in Figure 10L agree with the experimental data at [Na + ] i = [Na + ] o = 32 mM, [Ca 2+ ] i = 0, V o = 0, and V i = −20 mV (the intracellular membrane potential), where the partial Ca 2+ and Na + currents are denoted by the solid and dotted line, respectively. These two ionic currents show the anomalous mole fraction effect of the channel at nonequilibrium, i.e., trace concentrations of Ca 2+ ions effectively block the flow of abundant monovalent cations [212].
Free energies can be calculated by the electric and steric potentials [72] φ S2 = 1 4π 0 at the binding site S2 [215] on the atomic scale, where S2 also denotes Na + or K + (the site is occupied by a Na + or K + ), q j is the charge on the atom j in the protein given by PDB2PQR [219], p (r) = 1 + 77r/(27.7 + r) [119], r = |c j − c S2 |, c j is the center of atom j, A k is one of six symmetric surface points on the spherical S2, b = 3.6, and V S2 = 1.5v K + is a volume containing the ion at S2. We obtained ∆∆G = 5.26 kcal/mol [72] in accord with the MD result 5.3 kcal/mol [215], where G pore (Na + ) = 4.4, G bulk (Na + ) = −0.26, G pore (K + ) = −0.87, G bulk (K + ) = −0.27 kcal/mol, φ Na + = 7.5 k B T/e, v Na + v 0 S trc Na + = 0.23, φ K + = −1.93 k B T/e, v K + v 0 S trc K + = −0.59, and C B Na + = C B K + = 0.4 M. The crucial parameter in (60) is the ionic radius a S2 = 0.95 or 1.33 Å (also in |c j − A k |) that affects φ S2 very strongly but S trc S2 weakly. Another important parameter in (60) is the bulk void fraction Γ B that depends on the bulk concentrations of all ions and water and links the total energy of the ion at S2 to these bulk conditions measured very far away (∼10 6 Å) in the baths on the atomic scale.

Sodium Calcium Exchanger
The Na + /Ca 2+ exchanger (NCX) is the major cardiac mechanism that extrudes intracellular Ca 2+ across the cell membrane against its chemical gradient by using the downhill gradient of Na + [28].
The molecular basis of Na + /Ca 2+ interactions in NCX so striking to Lüttgau and Niedegerke [220] have been revealed by the cloning of NCX gene [221] and the structure of the ancient archaebacterial version NCX_Mj determined by Liao et al. [222]. Figure 12L illustrates NCX_Mj that consists of 10 transmembrane (TM) helices in which eight helices (TMs 2 to 5 and 7 to 10 labeled numerically in the figure) form a binding pocket of three putative Na + (green spheres) and one Ca 2+ (blue sphere) binding sites [222]. Structure of NCX_Mj consisting of ten transmembrane helices that form a binding pocket of three Na + (green spheres) and one Ca 2+ (blue sphere) binding sites [222]. Right (R): Schematic diagram of a cycle of Na + /Ca 2+ exchange in NCX consisting of five total potential states (TPS). Two Na + and one Ca 2+ ions enter the binding pocket in the outward-(TPS2 → TPS3 → TPS4) and inward-facing (TPS5 → TPS1) conformations, respectively. They exit in opposite conformations [70].
We developed a cyclic model of Na + /Ca 2+ exchange mechanism in NCX [70] using (60) to calculate five total (electric and steric) potential states (TPS) of various Na + and Ca 2+ ions shown in Figure 12R, where TPS1 and TPS4 are stable (with negative values) and TPS2, TPS3, and TPS5 are unstable (positive). Four extra sites in Figure 12R are determined empirically and close to entrance or exit locations of the binding pocket. The green and blue dots in the diagram represent Na + and Ca 2+ ions occupying the respective sites. Two Na + and one Ca 2+ ions enter the binding pocket in the outward-(TPS2 → TPS3 → TPS4) and inward-facing (TPS5 → TPS1) conformations, respectively. They exit in opposite conformations. The cycle consists of five steps.
Step 1: A conformational change is hypothetically activated [70] by a binding Ca 2+ at the blue site (S1) in TPS1 from inward-facing to outward-facing in TPS2.
Step 2: One Na + enters the binding pocket from the access site in TPS2 to the top Na + binding site (S2) in TPS3 followed by another Na + to the access site. These two coming Na + ions move the existing Na + ion from the middle Na + site (S3) to the bottom site (S4) by their Coulomb forces. TPS2 and TPS3 are unstable meaning that the two coming Na + ions have positive energies and are thus mobile. The selectivity ratio of Na + to Ca 2+ by NCX from the extracellular bath to the binding site S2 is C Na + (S2)/C Ca 2+ (S2) = 55.4 under the experimental conditions of the extracellular bath Na + o = 120 mM and Ca 2+ o = 1 µM [70].
Step 3: The vacant site S3 in TPS3 is a deep potential well with TP = −8.89 k B T/e that pulls the two unstable Na + ions to their sites in TPS4. Meanwhile, these two moving Na + and the stable Na + at S4 extrude the Ca 2+ (with unstable TP = 1.65) at S1 out of the pocket to become TPS4.
Step 4: Now, all three Na + ions in TPS4 are stable with negative TP and the vacant site S1 has an even deeper TP = −16.02. We conjecture that this TP value may trigger a conformational change from outward-facing in TPS4 to inward-facing in TPS5. The mechanism of conformational changes in NCX is yet to be studied.
Step 5: Furthermore, this large negative TP in TPS5 yields a remarkably large selectivity ratio of Ca 2+ to Na + by NCX from the intracellular bath to S1, i.e., C Ca 2+ (S1)/C Na + (S1) = 4986.1 at Ca 2+ i = 33 µM and Na + i = 60 mM. A coming Ca 2+ in TPS5 then extrudes two Na + ions out of the packet when it settles at S1 in stable TPS1.
Assuming that the total time T of an exchange cycle is equally shared by the 5 TPS, this model also infers that the stoichiometry of NCX is 3 5 T·2 Na + : 2 5 T·1 Ca 2+ = 3 Na + : 1 Ca 2+ in transporting Na + and Ca 2+ ions [70], which is the generally accepted stoichiometry (see reviews of Blaustein and Lederer [223] and Dipolo and Beaugé [224]) since the pivotal work of Reeves and Hale [225] and other subsequent experimental results. Please note that our model does not consider the electrogenic property of NCX [223], i.e., the driving force of the electric potential gradient.

Conclusions
We have covered a range of aspects of the fourth-order Poisson-Nernst-Planck-Bikerman theory from physical modeling, mathematical analysis, numerical implementation, to applications and verifications for aqueous electrolyte systems in chemistry and biology. The theory can describe many properties of ions and water in the system that classical theories fail to describe such as steric, correlation, polarization, variable permittivity, dehydration, mass conservation, charge/space competition, overscreening, selectivity, saturation, and more. All these properties are accounted for in a single framework with only two fundamental parameters, namely the dielectric constant of pure water and the correlation length of empirical choice. Ions and water have their physical volumes as those in molecular dynamic simulations. The theory applies to a system at both continuum and atomic scales due to the exact definition of the total volume of all ions, water molecules, and interstitial voids.
The most important features of PNPB are that (i) ions and water have unequal volumes with interstitial voids, (ii) their distributions are saturating of the Fermi type, (iii) these Fermi distributions approach Boltzmann distributions as the volumes tend to zero, and (iv) all the above physical properties appear self-consistently in a single model not separately by various models. Most existing modified Poisson-Boltzmann models consider ions of equal size and fail to yield Boltzmann distributions in limiting cases, i.e., the limit is divergent indicating that steric energies are poorly estimated. Numerous models for different properties such as steric, correlation, polarization, permittivity are proposed separately in the past.
We have shown how to solve 4PBik analytically and PNPB numerically. The generalized Debye-Hückel theory derived from the 4PBik model gives valuable insights into physical properties and leads to an electrolyte (analytical) equation of state that is useful to study thermodynamic activities of ion and water under wide ranges of composition, concentration, temperature, and pressure.
Numerically solving the fourth-order PNPB model in 3D for realistic problems is a challenging task. There are many pitfalls that one must carefully avoid in coding. For that reason, we have particularly mentioned some methods for handling the convergence issues of the highly nonlinear PNPB system of partial differential equations and the discretization problems concerning the complicate interface between molecular and solvent domains and the Scharfetter-Gummel stability condition to ensure positivity of numerical concentrations and current preservation.
Finally, we have shown novel results obtained by PNPB for chemical and biological systems on ion activities, electric double layers, gramicidin A channel, L-type calcium channel, potassium channel, and sodium calcium exchanger. These results agree with experiments or molecular dynamics data and show not only continuum but also atomic properties of the system under far-field conditions. The fourth-order PNPB model is consistent and applicable to a great variety of systems on a vast scale from meter to Angstrom.

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

Abbreviations
The following abbreviations are used in this manuscript: