Poisson-Fermi Formulation of Nonlocal Electrostatics in Electrolyte Solutions

We present a nonlocal electrostatic formulation of nonuniform ions and water molecules with interstitial voids that uses a Fermi-like distribution to account for steric and correlation effects in electrolyte solutions. The formulation is based on the volume exclusion of hard spheres leading to a steric potential and Maxwell's displacement field with Yukawa-type interactions resulting in a nonlocal electric potential. The classical Poisson-Boltzmann model fails to describe steric and correlation effects important in a variety of chemical and biological systems, especially in high field or large concentration conditions found in and near binding sites, ion channels, and electrodes. Steric effects and correlations are apparent when we compare nonlocal Poisson-Fermi results to Poisson-Boltzmann calculations in electric double layer and to experimental measurements on the selectivity of potassium channels for K+ over Na+. The present theory links atomic scale descriptions of the crystallized KcsA channel with macroscopic bulk conditions. Atomic structures and macroscopic conditions determine complex functions of great importance in biology, nanotechnology, and electrochemistry.

ical properties of electrolyte solutions in a wide range of applications in electrochemistry, biophysics, colloid science, and nanofluidics [1][2][3][4][5][6]. For over a century, a great deal of effort has been devoted to improving the Poisson-Boltzmann (PB) theory of continuum models for a proper description of steric (or finite size) and correlation (or nonlocal screening, polarization) effects in electrolytes [7,8]. We present a continuum model with Fermi-like distributions and global electrostatic screening of nonuniform ions and water molecules to describe the steric and correlation effects, respectively, in aqueous electrolyte solutions.
For an electrolytic system with K species of ions, the entropy model proposed in [9] treats all ions and water of any diameter as nonuniform hard spheres and regards the water as the (K + 1) th species. It then includes the voids between these hard spheres as the (K + 2) th species so that the total volume V of the system can be calculated exactly by the identity where V K+2 denotes the total volume of all the voids, v i = 4πa 3 i /3 that gives the volume of each sphere with radius a i , and N i is the total number of the i th species particles. In the bulk solution, we have the bulk concentrations C B i = N i V and the bulk volume fraction of voids Γ B = V K+2 V . Dividing the volume identity (1) by V , Γ B = 1 − K+1 i=1 v i C B i is expressed in terms of nonuniform v i and C B i for all particle species. If the system is spatially inhomogeneous with variable electric or steric fields, as in most biological and technological systems, the bulk concentrations then change to concentration functions C i (r) that vary with positions, and differ from their constant values C B i at location r in the solvent domain Ω s . Consequently, the void volume fraction becomes a function Γ(r) For an electrolyte in contact with electrodes or containing a charged protein, an electric field E(r) in the solvent domain Ω s is generated by the electrodes, ionic (free) charges with a displacement field D(r), and bound charges of polar water with a polarization field P(r).
In Maxwell's theory, these fields form a constitutive relation that yields the Maxwell's equation ∇ · D(r) = ρ(r) = K i=1 q i C i (r), ∀r ∈ Ω s , where ǫ 0 is the vacuum permittivity and q i is the charge on each i species ion [10]. The electric field E(r) is thus screened by water (in what might be called Bjerrum screening) and ions (in Debye screening) in a correlated manner that is usually characterized by a correlation (screening) length λ [4,6,11]. 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 screening kernel G(r − r ′ ) = e −|r−r ′ |/λ 4π|r−r ′ | [1], which is called Yukawa-type kernel in [4,12], and satisfies the partial differential equation (PDE) [12] − ∆G(r − r ′ ) + 1 in the whole space R 3 , where ∆ = ∇ · ∇ is the Laplace operator with respect to r and thus describes a local potential of free ions according to the Poisson equation where ǫ s is a dielectric constant in the solvent domain. We introduce a global electric potential φ(r) of the screened electric field E(r) as a convolution of the local potential φ(r ′ ) with the global screening kernel G(r − r ′ ) in the expression However, it would be too expensive to calculate φ(r) using this equation. Multiplying Eq.
If λ = 0, we recover the standard Poisson equation (4) and the classical field relation P = ǫ 0 (ǫ s − 1)E with the electric susceptibility ǫ s − 1 (and thus the dielectric constant ǫ s ) if water is treated as an isotropic and linear dielectric medium [10].
We introduce a Gibbs free energy functional for the system as an average volume, and L −1 is the inverse of the self-adjoint positive linear operator L = ǫ s λ 2 ∆∆ − ǫ s ∆ [12]. Taking the vari- for all i = 1, · · · , K + 1 (ions and water), where Boltzmann constant, and T is an absolute temperature. The distribution is of Fermi type since it saturates. All concentration functions C i (r) < 1 v i [9], i.e., C i (r) cannot exceed the maximum value 1/v i for any arbitrary (or even infinite) potential φ(r) in the domain Ω s .
In these Fermi distributions, it is impossible for a particle 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 yield S trc (r) = −∞ and hence v i C i (r) = 0, a contradiction. For this reason, we must include the void as a separate species if water and ions are all treated as hard spheres [9]. Here we do represent water and ions as spheres. Our approach allows other shapes to be used as well.
(i) The Fermi-like distribution of uniform spherical ions with voids in ionic liquids was first derived by Bazant et al. [6,7] using Bikerman's lattice model [13]. The entropy functional in [6] involves a reciprocal of a minimum volume v with a volume fraction Φ that is an empirical fitting parameter and may have to be unrealistically large to fit experimental data in some applications [7]. It is shown in [9] that the entropy functional in [6] does not directly yield classical Boltzmann distributions C i (r) = C B i exp (−β i φ(r)) as v → 0. It can be easily seen from (9) that the entropy functional F en (C) in Eq. (8) consistently yields Boltzmann distributions as v i → 0 for all i. Our derivation of F en (C) does not employ any lattice models but simply uses the volume equation (1). The functional F en (C) is a new modification of that in [9], where the classical Gibbs entropy is now generalized to include all species (ions, water, and voids) in electrolytes with the same entropy form. In fact, our Γ(r) is an analytical extension of the void fraction 1 − Φ in Bikerman's excess chemical potential [7], where all volume parameters v i (including the bulk fraction Γ B ) are physical not empirical. The adjustable parameter in our model is the correlation length λ ≈ 2a i depending on the ionic species i of interest [9,14]. The PF model was first proposed in [7] without derivation and has been shown to produce the results that are not only comparable to molecular dynamics (MD) simulation or experimental data but also provide insight into nonlinear properties of concentrated electrolytes and ionic liquids [15]. Here, we formally derive the PF model for general electrolytes using a hard-sphere instead of lattice model with the steric potential S trc (r) first introduced in [14]. As compared with lattice models in [7] and demonstrated, for example, in [9,16], hard-sphere models significantly improve the agreement between simulation and experiment.
However, the physical origin of the PDE is different. In [11], the global convolution is performed only on the charge density of point-like counter ions in contrast to the potential φ(r) by Eq. (4) that is generated by all spherical ions. In [6], a more general derivation for electrolytes or ionic liquids with steric effects is given from a free energy function of a gradient expansion of nonlocal electrostatic energy in terms of ∆ φ. Eq. (7) corresponds to the first term in the expansion. Here, the fourth-order PDE is derived directly from Maxwell's equation with the Yukawa screening kernel. Our result does not depend on the convergence properties of an expansion of nonlocal electrostatic energy.
(iii) Eq. (7) defines an dielectric operator ǫ = ǫ s ǫ 0 (1 − λ 2 ∆) that in turn implicitly yields a dielectric function ǫ(r) as an output of solving Eq. (7) [6,9]. A physical interpretation of the operator was first introduced in [6] to describe the nonlocal permittivity in a correlated ionic liquid. The exact value of ǫ(r) at any r ∈ Ω s cannot be obtained from Eq. (7) but can be approximated by the simple formula ǫ(r) ≈ ǫ 0 + C H 2 O (r)(ǫ s − 1)ǫ 0 /C B H 2 O since the water density function C H 2 O (r) = C K+1 (r) is an output of Eq. (9). This formula is only for visualizing (approximately) the profile of ǫ or ǫ. It is not an input of calculation. The input is the operator ǫ = ǫ s ǫ 0 (1 − λ 2 ∆) (or the correlation length λ).
(iv) The factor v i /v 0 multiplying S trc (r) in Eq. (9) is a modification of the unity used in our previous work [9]. The steric energy − v i v 0 S trc (r)k B T [9] of a type i particle depends not only on the emptiness (Γ(r) = 1 − K+1 i=1 v i C i (r)) (or equivalently crowding) at r but also on the volume v i of each type of particle. If all v i are equal (and thus v i = v 0 ), then all particle species at any location r ∈ Ω s have the same steric energy and the uniform particles are indistinguishable in steric energy. The steric potential is a mean-field approximation of Lennard-Jones (L-J) potentials that describe local variations of L-J distances (and thus empty voids) between every pair of particles. L-J potentials are highly oscillatory and extremely expensive and unstable to compute numerically.
(v) The global convolution in Eq. (5) may seem similar to those in [4,12] but it is not.
The Poisson equation (4) takes the place of the Fourier-Lorentzian (FL) equation -an integro-differential equation -in the previous work [4,12] in which the local potential φ(r ′ ) in Eq. (5) is replaced by a global electric potential Φ(r ′ ). Moreover, the factor 1/λ 2 in Eq. A in contrast to an edge length of 7.5Å of cubical ions in [17], e is the proton charge, and ǫ s = 80. The multivalent ions A 4− represent large polyanions adsorbed onto a charged Langmuir monolayer in experiments [17]. The dotted curve in Fig. 1 is similar to that of mPB in [17] and was obtained by the PF model with the size effect but without voids and correlations, i.e., λ = 0, V K+2 = 0 (no voids), and v K+1 = v H 2 O = 0 (water is volumeless and We next show the size effect of different ions with voids using the crystal structure of the potassium channel KcsA (PDB [20] ID 3F5W [21]) as shown in Fig. 3, where the spheres denote five specific cation binding sites (S0 to S4) [22]. The crystal structure with a total of N = 31268 charged atoms is embedded in the protein domain Ω p while the binding sites are in the solvent domain Ω s . The exquisite selectivity of K + (with a K + = 1.33Å) over other cations such as Na + (a Na + = 0.95Å) by potassium channels is an intriguing quest in channology. It can be quantified by the free energy (G) difference of K + and Na + in the pore and in the bulk solution [22,23]. Experimental measurements [23] showed that the relative free energy [22] ∆∆G(K + → Na + ) = G pore (Na is greater than zero in the order of 5-6 kcal/mol unfavorable for Na + . The electric and steric potentials at the binding site S2 (as considered in [22]) can be calculated on the atomic scale using the following algebraic formulas where S2 = Na + or K + (the site is occupied by a Na + or a K + ), q j is the charge on the atom j in the protein given by PDB2PQR [24], ǫ p (r) = 1 + 77r/(27.7 + r) [25], 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. The crucial parameter in (12) 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. We obtained ∆∆G = 5.26 kcal/mol in accord with the MD result in [22], where G pore (Na + ) = 4.4, G bulk (Na + ) = −0.26 [26], G pore (K + ) = −0.87,