Analysis of a Poisson-Nernst-Planck-Fermi system for charge transport in ion channels

A modified Poisson-Nernst-Planck system in a bounded domain with mixed Dirichlet-Neumann boundary conditions is analyzed. It describes the concentrations of ions immersed in a polar solvent and the correlated electric potential due to the ion--solvent interaction. The concentrations solve cross-diffusion equations, which are thermodynamically consistent. The considered mixture is saturated, meaning that the sum of the ion and solvent concentrations is constant. The correlated electric potential depends nonlocally on the electric potential and solves the fourth-order Poisson-Fermi equation. The existence of global bounded weak solutions is proved by using the boundedness-by-entropy method. The novelty of the paper is the proof of the weak--strong uniqueness property. In contrast to the existence proof, we include the solvent concentration in the cross-diffusion system, leading to a diffusion matrix with nontrivial kernel. Then the proof is based on the relative entropy method for the extended cross-diffusion system and the positive definiteness of a related diffusion matrix on a subspace.


Introduction
The modeling of the transport of ions through biological channels is of fundamental importance in cell biology. Several strategies have been developed in past decades, using molecular or Brownian dynamics or the Poisson-Nernst-Planck theory. This theory relies on the assumptions that the dynamics of ion transport is based on diffusion and electrostatic interaction only and that the solution is dilute. However, the presence of narrow channel pores requires a more sophisticated modeling. In particular, the ion size is not small compared to the biological channel diameter, and many-particle interactions due to the confined geometry need to be taken into account. In this paper, we analyze a modified Poisson-Nernst-Planck system modeling ion-water interactions and finite ion size constraints. We prove the existence of global weak solutions and, as the main novelty, the weak-strong uniqueness property using entropy methods.
The unknowns are the ion concentrations (or volume fractions) u i (x, t) of the ith ion species and the correlated electric potential Φ(x, t). The solvent concentration (or volume fraction) u 0 (x, t) is given by u 0 = 1 − n i=1 u i , which means that the mixture is saturated. Equations (1) are cross-diffusion equations with the fluxes J i and the reaction rates r i (u). The parameters are the diffusivities D i > 0 and the valences z i ∈ Z. Equation (2) is the Poisson-Fermi equation with the scaled Debye length λ > 0, the correlation length ℓ > 0, and the given background charge density f (x). We assume that the domain is isolated on the Neumann boundary, while the concentrations and the electric potential are prescribed on the Dirichlet boundary.
In the following, we discuss definition (1) of the fluxes and equation (2) for the correlated electric potential.
We recover the classical Poisson-Nernst-Planck equations if u 0 = const. and ℓ = 0. In this situation, we can write J id i = −D i u i ∇µ id i with the electrochemical potential µ id i = log u i + z i Φ of an ideal dilute solution. In concentrated solutions, the finite size of the ions needs to be taken into account, expressed by the excess chemical potential µ ex i , so that the electrochemical potential becomes µ i = µ id i + µ ex i . Bikerman [4] suggested the choice µ ex i = − log(1 − n i=1 u i ) = − log u 0 ; also see [1,Sec. 3.1.2]. Then J i = −D i u i ∇µ i coincides with the flux adopted in our model (1). Note that solving µ i = log(u i /u 0 ) + z i Φ for the concentrations, we find that the ion profiles obey the Fermi-Dirac statistics , i = 1, . . . , n.
Then, given µ i and Φ, the bounds 0 ≤ u i ≤ 1 are automatically satisfied. Other choices of the excess chemical potential were suggested in [3,Sec. 2.1].
In the literature, there exist also other approaches to define the fluxes J i under finite size constraints. The diffusion limit of an on-lattice model, which takes into account that neighboring sites may be occupied (modeling size exclusion), was performed in [6], analyzed where the dielectricity ε 0 = λ 2 , instead of being constant, depends on u. This assumption is based on the experimental observation that the dielectric response of water decreases as ion concentrations increase [8]. Thus, ∂ε 0 /∂u i < 0, showing that the ion-water interaction energy is always nonnegative. Finite ion size effects are modeled in [20] by including an approximation of the Lennard-Jones potential in the energy functional, leading to Assuming that (a ij ) is positive definite, the global existence of weak solutions for two species was proved in [18]. For the analysis of the stationary equations, see [14]. Finally, excluded volume effects can be included by considering nonlinear diffusivities D i (u i ) = 1 + αu i , where α > 0 is a measure of the volume exclusion interactions [5]. Our model has the advantage of being consistent with the thermodynamical model [11] assuming that the diffusion matrix is diagonal, D ij = D i δ ij , and that the solvent is neutral, z 0 = 0. The interaction of the ions with polar solvents like water is modeled by the potential in (2). Indeed, let φ be the electric potential of free ions, given by −λ 2 ∆φ = ρ, where ρ is the total charge density. Then the correlated potential Φ = ℓ −2 Y ℓ * φ is the convolution between the Yukawa potential Y ℓ (x) = (|x|/ℓ) −1 exp(−|x|/ℓ) [21] and the electric potential, where ℓ > 0 is the correlation length of the screening by ions and water [24]. As this potential satisfies −ℓ 2 ∆Φ + Φ = φ, we recover (2) with ρ = n j=1 z j u j + f (x). Thus, the Poisson-Fermi equation (2) includes finite ion size effects and polarization correlations among water molecules. It generalizes the fourth-order differential permittivity operator of [22] and the nonlocal permittivity in ionic liquids of [2]. If there are no correlation and polarization effects (ℓ = 0), we recover the standard Poisson equation for the electric potential. The expression ε 0 = λ 2 (ℓ 2 ∆ − 1) can be interpreted as a dielectric differential operator.
1.2. Entropy structure. System (1) can be written as a cross-diffusion system with a diffusion matrix which is neither symmetric nor positive definite. This issue is overcome by exploiting the entropy (or free energy) structure and using the boundedness-by-entropy method [19]. The free energy associated to (1)-(2) is given by [2] The energy density h(u) consists of the internal, free-ion electric, and correlation electric energies. The free energy allows us to formulate equations (1) as a diffusion system with a positive semidefinite diffusion matrix. Indeed, we introduce the electrochemical potentials where ∂h/∂u i denotes the variational derivative of h with respect to u i (see [15,Lemma 7]) and u D 0 = 1 − n i=1 u D i . As in [15], we split the electrochemical potentials into the entropy variables w i and the boundary contributions w D i by Then equations (1) can be written as is interpreted as a function of w = (w 1 , . . . , w n ) and Φ according to .
The advantage of formulation (9) is that the new diffusion matrix B = (B ij ) is symmetric and positive semidefinite. Observe that system (9) is of degenerate type since u i = 0 is possible, and det B = 0 in this case. The formulation in terms of entropy variables has the further advantage that the ion concentrations u i , defined by (10), are nonnegative and satisfy n i=1 u i ≤ 1, thus fulfilling the saturation assumption.
1.3. Main results. We introduce the simplex D = {u = (u 1 , . . . , u n ) ∈ (0, 1) n : n i=1 u i < 1} and set Ω T = Ω × (0, T ). The following hypotheses are imposed: (H5) Reaction rates: r i ∈ C 0 ([0, ∞) n ; R) for i = 1, . . . , n, and there exists C r > 0 such that for all u ∈ L ∞ (Ω T ; D) and Φ, given by (2) and (5), The restriction to three space dimensions in Hypothesis (H1) is not needed. It can be removed by regularizing the Poisson-Fermi equation (2) to ensure that Φ ∈ L ∞ (Ω); see Remark 3 below. In Hypothesis (H4), it is sufficient to define the boundary data on Γ D . We have extended them to Ω with the special extention of Φ D , fulfilling the fourth-order elliptic problem (11). This extension is needed to be consistent with the definition of the free energy and the entropy variables; see [15,Lemma 7]. The bound in Hypothesis (H5) is needed to derive gradient bounds on the concentrations from the free energy inequality; see (16) below. Since ∂h/∂u i contains the logarithm, r i (u) needs to cancel the singularity in ∂h/∂u i at u i = 0. The contribution of the potential Φ − Φ D in ∂h/∂u i is treated by applying the Poincaré inequality and observing that Ω |∇(Φ − Φ D )| 2 dx ≤ CH(u) for some C > 0. Therefore, we need the integrated version (12) instead of the pointwise inequality assumed in [19,Sec. 1.4].
We introduce the test spaces Our first main result is as follows.
; L 2 (Ω)), Φ ∈ L 2 (0, T ; H 2 (Ω)), log u 0 ∈ L 2 (0, T ; H 1 (Ω)), for all φ i ∈ L 2 (0, T ; H 1 D (Ω)) and θ ∈ L 2 (0, T ; H 2 D,N (Ω)), where J i is given by (1) and ·, · is the dual product between in Ω), the solution satisfies for 0 < s < t < T the free energy inequality The energy dissipation is understood in the sense The assumption of thermal equilibrium at the Dirichlet boundary, also required in [15], is needed to avoid expressions involving ∇w D i in the free energy inequality. Thus, this condition, together with vanishing reactions, is natural to obtain the monotonicity of the free energy. We refer to Remark 9 for the case with reaction rates.
The proof of Theorem 1 is, similarly as in [15], based on an approximation procedure, where we regularize (9) by an implicit Euler approximation and higher-order terms in the entropy variables. The uniform estimates that are needed to perform the de-regularization limit are derived from the free energy inequality, which (without regularization) reads as recalling definition (8) of w i , and we can conclude by Gronwall's lemma. The free energy dissipation term on the left-hand side can be estimated from above by (see Lemma 6) The last term is bounded by the electric energy part in H(u), thus giving H 1 (Ω) bounds for u i for i = 0, . . . , n and log u 0 . Compared to [15], we obtain gradient estimates for all the ion concentrations, but we have to deal with the singular term ∇ log u 0 in (1). Moreover, compared to [13], where a similar Nernst-Planck system (with ℓ = 0) was investigated, we do not need any positivity condition on the initial solvent concentration.
While the existence proof relies on standard entropy methods, we need a new idea to prove the weak-strong uniqueness result. The uniqueness of weak solutions is an intricate problem. A uniqueness result for (1) with the fluxes (6) was shown in [15] for the case D i = D and z i = z for all i. In this simplified situation, the solvent concentration solves a Poisson-Nernst-Planck system for which the uniqueness of bounded weak solutions can be proved by a combination of L 2 (Ω) estimates and Gajewski's entropy method. This strategy cannot be used for our system. In fact, we need the H −1 (Ω) method and a strong regularity condition for ∇Φ, which restricts the geometry of the Dirichlet-Neumann boundary conditions; see Remark 14. Therefore, we do not aim to prove the uniqueness of weak solutions but the weak-strong uniqueness property only, which has the advantage that we may allow for different coefficients D i and z i . The weak-strong uniqueness property means that any weak solution to system (1)-(5) coincides with a strong solution emanating from the same initial conditions as long as the latter exists. We say that (ū,Φ) is a strong solution to (1)-(5) if it is a weak solution and Our second main result is contained in the following theorem.
Theorem 2 (Weak-strong uniqueness). Let the Dirichlet boundary data be in thermal equilibrium in the sense of Theorem 1 and let r i = 0 for i = 1, . . . , n. Let (u, Φ) be a weak solution and (ū,Φ) be a strong solution to (1)- (5).
If the reaction rates are Lipschitz continuous and satisfy some sign conditions, Theorem 2 still holds; see Remark 15. The condition that the Dirichlet boundary data is in thermal equilibrium is actually not needed, since in contrast to (15), the terms involving ∇w D i cancel out in the computations for the relative free energy which can be identified as the Bregman distance of the free energy. The key idea of the proof of Theorem 2 is to consider the solvent concentration u 0 as an independent variable and to formulate the parabolic equations for the extended concentration vector U = (u 0 , u 1 , . . . , u n ), leading to This situation is similar to the Maxwell-Stefan system; see [17]. The time derivative of the relative free energy equals and K 2 contains differences like U i −Ū i and Φ −Φ. The properties of the matrices (A ij ) and (Q ij ) imply that where P L is the projection on L and Y i = √ u i ∇ log(u i /ū i ), as well as for any δ > 0, Consequently, choosing δ > 0 sufficiently small, for some constant C > 0. Since the initial data of u andū coincide, we have H((u, Φ)(t)| (ū,Φ)(t)) = 0 and finally u(t) =ū(t) and Φ(t) =Φ(t) for all t > 0. The idea to consider the parabolic system for the extended solution vector U = (u 0 , . . . , u n ) instead of u = (u 1 , . . . , u n ) is the main novelty of this paper. The article is organized as follows. The proof of Theorem 1 is presented in Section 2, while Section 3 contains the proof of Theorem 2. We make some remarks on the uniqueness of solutions in Section 4.

Proof of Theorem 1
We assume throughout this section that Hypotheses (H1)-(H5) hold.
2.1. Solution of an approximate system. We define the approximate problem by the implicit Euler scheme and using a higher-order regularization. Let T > 0, N ∈ N, τ = T /N, and m ∈ N with m > d/2. We assume that u D i ≥ η > 0 for i = 0, . . . , n. Then Since the entropy variables are not needed in the weak formulation (13)- (14), we can pass to the limit η → 0 at the end of the proof, thus requiring only is the unique solution to (2) with u 0 j instead of u j on the right-hand side and satisfying the corresponding boundary conditions in (4)- (5). We wish to find a solution v k ∈ X := H m (Ω; . . , r n (u)), and Thanks to the higher-order regularization, we obtain approximate solutions w k : Remark 3. Adding a higher-order regularization to the Poisson-Fermi equation (18), we may obtain Φ k ∈ L ∞ (Ω) by a Sobolev embedding similarly as for w k . This allows us to remove the restriction d ≤ 3 in Hypothesis (H1).

Lemma 4. There exists a unique solution
The proof is similar to that one of Lemma 5 in [15], therefore we give a sketch only. Let y ∈ L ∞ (Ω; R n ) and σ ∈ [0, 1]. Let Φ k ∈ H 2 (Ω) be the unique solution to in Ω subject to the boundary conditions (5). This follows from the fact that the function (x, Φ) → u i (w(x), Φ) is bounded with values in (0, 1) and Lipschitz continuous in Φ. By the Lax-Milgram lemma, there exists a unique solution v ∈ X to the linear problem Indeed, as B is positive semidefinite, the left-hand side is coercive in H m (Ω; R n ). This defines the fixed-point operator S : Then S(y, 0) = 0, S is continuous and, because of the compact embedding H m (Ω; R n ) ֒→ L ∞ (Ω; R n ), also compact. Using φ = v as a test function in (19), standard estimates lead Hence, all fixed points of S(·, δ) are uniformly bounded in L ∞ (Ω; R n ). We infer from the Leray- 2.2. Uniform estimates. We deduce estimates uniform in (ε, τ ) from the following free energy inequality.
where H is defined in (7) and C r > 0 is introduced in Hypothesis (H5). (17). Using the generalized Poincaré inequality to estimate the ε-regularization and Hypothesis (H5) to estimate the reaction rates, we find that It follows from the convexity of the function g(u) = n i=0 u i u D i log(s/u D i )ds and the Poisson-Fermi equation (2) as in [15,Section 2] that Inserting the definition B ij (w k , Φ k ) = D i u k i δ ij , we infer from Young's inequality that Collecting these estimates and observing that u k i ≤ 1 concludes the proof. We sum (20) over k = 1, . . . , j, and, assuming τ < 1/C r , apply the discrete Gronwall inequality [9]: where C(T ) > 0 does not depend on (ε, τ ). We still need to bound the second term on the left-hand side from below.
where C > 0 depends on (D i ) and (z i ).
Proof. We infer from Young's inequality and the bound u k i ≤ 1 that The first term on the right-hand side is rewritten as using u k 0 ≤ 1 in the last step. Since the free energy is bounded from below, we conclude the following uniform bounds. Lemma 7. There exists C > 0 not depending on (ε, τ ) such that for i = 1, . . . , n, Proof. The inequality The H 2 (Ω) bound for Φ k follows immediately from the Poisson-Fermi equation as its right-hand side is bounded in L 2 (Ω). The H 1 (Ω) bound for log u k 0 is a consequence of the L 2 (Ω) bound for ∇ log u k 0 and the Poincaré inequality, using the fact that log u D 0 ∈ L 2 (Ω) by Hypothesis (H4).
√ ε w where i = 1, . . . , n. We also need a uniform bound for the discrete time derivative. Lemma 8. There exists a constant C > 0 independent of (ε, τ ) such that for all i = 1, . . . , n, By a density argument, this inequality holds for all φ i ∈ L 2 (0, T ; X), showing the desired bound for the discrete time derivative of u We claim that ∇ log u (τ ) 0 ⇀ ∇ log u 0 weakly in L 2 (Ω T ). It follows from (24) that (for a subsequence) ∇ log u (τ ) 0 ⇀ v weakly in L 2 (Ω T ). We need to identify v = ∇ log u 0 . We know that (again for a subsequence) u shows that log u (τ ) 0 → log u 0 strongly in L 2 (Ω T ). Hence, we conclude that v = ∇ log u 0 , proving the claim.
Remark 9. Let the reaction rates r i : D → R be Lipschitz continuous and quasi-positive, i.e. r i (u) ≥ 0 for all u ∈ D with u i = 0. We assume that the total reaction rate is nonnegative, i.e. n i=1 r i (u) ≤ 0 for all u ∈ D, and that r i (u) log u i = 0 if u i = 0. We claim that the free energy inequality becomes This inequality follows from (26) after including the reaction rates and taking the limit Indeed, the strong limit u in Ω T if u i > 0. If u i = 0, by assumption, we have r i (u) log u i = 0 and therefore r i (u (τ ) ) log u (τ ) i → r i (u) log u i a.e. in Ω T as well. Moreover, r i (u) log u i is bounded. Hence, by dominated convergence, r i (u (τ ) ) log u (τ ) i → r i (u) log u i strongly in L 1 (Ω T ), and the claim follows.

Proof of Theorem 2
Let (u, Φ) be a weak solution and (ū,Φ) be a strong solution to (1)- (5). In this section, we interpret H(u) and H(ū) as functionals depending on u = (u 0 , . . . , u n ) andū = (ū 0 , . . . ,ū n ). This notation is only needed to determine the variational derivative of H and will not lead to any confusion in the following computations. We split the lengthy proof in several steps.
Step 1: Calculation of the time derivative of H(u, Φ|ū,Φ). In the following, we write is the variational derivative of H 1 atū in the direction of u −ū (similarly for H ′ 2 (Φ)(Φ −Φ)). We compute the time derivative of H 1 (u|ū), split the sum over i = 0, . . . , n into i = 0 and the sum over i = 1, . . . , n, and Next, we insert equation (1) for u i andū i and use (∂h A similar computation for H 2 (Φ|Φ) leads to where we abbreviatedw i = log(ū i /ū 0 ) + z iΦ . As u i is only nonnegative, the expression ∇ log u i may be not integrable. Therefore, we define ∇ log(u i /u 0 ) := (2∇ √ u i − √ u i ∇ log u 0 )/ √ u i if u i > 0 and ∇ log(u i /u 0 ) := 0 else. This expression may be not In a similar way, we define ∇ log(u i /ū i ), which is possible sinceū i is strictly positive, and we have We insert the free energy inequality (15), namely and rearrange the terms, At this point, we observe that the terms involving ∇w D i cancel even if ∇w D i does not vanish, since they also appear in the free energy inequality (15).
The treatment of I 1 and I 2 is more delicate.
Step 3: Estimation of I 1 . We write I 1 = I 11 + I 12 + I 13 , where It follows from 0 ≤ u i ≤ 1 that |Q ij | ≤ C and consequently The matrix A is not positive definite since u i = 0 is possible. We show that the matrix G, defined by . The corresponding positive bound will help us to estimate I 13 and I 2 .
Proof. Let G ij = A ij / √ u i u j if u i u j > 0 and G ij = 0 else. Recall that by definition, In this case, If u i = 0 or u j = 0, either Y i = 0 or Y j = 0 and hence, the previous expression vanishes. Therefore, We know that u 0 > 0 a.e. in Ω T . A straightforward computation shows that ran G = L, implying that ker G = L ⊥ . We claim that the matrix G satisfies for all Y ∈ R n+1 , Indeed, introduce first the matrix Let Y ∈ R n+1 and z i : This implies that is integrable in Ω T since ∇ log u 0 ∈ L 2 (Ω T ), and √ u j Y j ∈ L 2 (Ω T ). Therefore, we can integrate inequality (33), and in view of GY = G(P L Y ) (since ker G = L ⊥ ), this shows that which finishes the proof.
For any ε > 0, there exists C(ε) > 0 such that Proof. We take into account the structures of the matrices A and Q: This gives Next, we calculate for i = 0, . . . , n, where we used n j=0 ∇u j = 0 and n j=0ū j ∇ logū j = 0. Hence, We split Y i = (P L Y ) i + (P L ⊥ Y ) i in I 13 , which leads to An application of Young's lemma finishes the proof.
The previous lemmas show that Step 4: Estimation of I 2 . We split I 2 = I 21 + I 22 , where Lemma 12. For any ε > 0, there exists C(ε) > 0 such that Proof. Recalling that Y i = √ u i ∇ log(u i /ū i ), we reformulate I 21 as All rows of the matrix (A ij /u i −Ā ij /ū i ) vanish except the first one, This shows that Since ∇ logū i is bounded in L ∞ (Ω T ), we can bound the first term in I 21 : where ε > 0 is arbitrary. To estimate the second term in I 21 , we use (34) and the elementary inequality ( n i=0 |u i −ū i |) 2 ≤ (n + 1) n i=0 |u i −ū i | 2 : The lemma follows after inserting (38) and (39) into (37).
Proof. All entries of the matrix (Q ij /u i −Q ij /ū i ) vanish except the element Q 00 /u 0 − Q 00 /ū 0 = − n i=1 D i z i (u i /u 0 −ū i /ū 0 ). This leads to It follows from (34) that Hence, Young's inequality completes the proof.
Remark 15 (Weak-strong uniqueness in the presence of reaction terms). We claim that Theorem 2 holds for reaction rates r i : D → R, which are Lipschitz continuous and quasipositive (i.e. r i (u) ≥ 0 for all u ∈ D with u i = 0) such that the total reaction rate is nonnegative, i.e. n i=1 r i (u) ≤ 0 for all u ∈ D, and that r i (u) log u i = 0 if u i = 0. Proceeding as in Step 1 of the proof of Theorem 2 and taking into account Remark 9, we need to estimate additionally the expression The assumptions on r i imply that r i (u) log u i is integrable. Therefore, following [12, p. 202f], We deduce from 0 ≥ log z − z + 1 ≥ −|z − 1| 2 / min{1, z} for z > 0 that where we used in the last step the assumption n i=1 r i (u) ≤ 0. Furthermore, by the Lipschitz continuity of r i , the Poincaré inequality, and the elliptic estimate for the Poisson-Fermi equation, Thus, estimate (41) is still valid with another constant, and Theorem 2 follows.