Existence of Approximate Solutions for Modified Poisson Nernst-Planck Describing Ion Flow in Cell Membranes

Abidha Monica Gwecho¹, Shu Wang², Onyango Thomas Mboya³

¹College of Applied Sciences, Beijing University of Technology, Beijing, China.
²School of Mathematics, Beijing University of Technology, Beijing, China.
³School of Mathematics and Actuarial sciences, Technical University of Kenya, Nairobi, Kenya.

DOI: 10.4236/ajcm.2020.103027                

Abstract

Dynamics of ions in biological ion channels has been classically analyzed using several types of Poisson-Nernst Planck (PNP) equations. However, due to complex interaction between individual ions and ions with the channel walls, minimal incorporation of these interaction factors in the models to describe the flow phenomena accurately has been done. In this paper, we aim at formulating a modified PNP equation which constitutes finite size effects to capture ions interactions in the channel using Lennard Jonnes (LJ) potential theory. Particularly, the study examines existence and uniqueness of the approximate analytical solutions of the mPNP equations, First, by obtaining the priori energy estimate and providing solution bounds, and finally constructing the approximate solutions and establishing its convergence in a finite dimensional subspace in L², the approximate solution of the linearized mPNP equations was found to converge to the analytical solution, hence proof of existence.

Keywords

Lennard Jonnes Potential, Finite Size Effects, Ion Channel, Modified PNP

1. Introduction

Biological cells are composed of proteins arranged in folded chains of amino acids to form ionic channels that are nanoscopic water-filled pores to perform the role of controlling transport of ions in cell membranes. The channel maintains correct ion composition and balance in cells that is crucial in their survival and numerous functions of propagating life such as, energy conversion, drug delivery, secretion, among others. These functions are varied and enabled by ability of the cells to carry strong and steeply varying distribution of permanent charges depending on combination of the nanotubes and prevalent physiological conditions.

Ion channels are characterized by their functioning, some are known to exhibit complex switching properties similar to electronic devices, others have the ability to selectively transport or block a particular ion species, while others have no selectivity, see [1]. To develop deeper understanding of the processes in the channel both analytical and empirical investigation are critical. Numerical experiments approximate transport through biological channels to determine amongst other things structures and conductance of ion channels as a means of minimizing cost and complementation of empirical findings, see [2]. Poisson Nernst Planck (PNP) equations for a long period of time has been adopted as a classical mathematical model and analysis tool of choice for studying ion flow. The model couples electrostatics with diffusion process as a popular theoretical method to robustly simulate ion channel systems. However, the major drawback of PNP model is that it neglects finite size effects in biological channel systems resulting into significant inaccuracies.

Incorporating electrostatic interaction of ions and finite size effects particularly in narrow regions has been suggested, investigated and determined to impact in reduction of error in solution approximation. Energy Variational Approach was used by [3] to derive an accurate generalized Poisson Nernst Planck-Navier Stokes (PNP-NS) system which characterizes interactions of charged fluid and mutual friction between the crowded charged particles. A general method was thereafter developed to show that the system is globally asymptotically stable under small perturbation around a constant equilibrium. Subsequently, [4] derived mPNP system which includes an extra dissipation due to effective velocity differences between ion species, then using Galerkins method and Schauders fixed point theorem local existence theorem of the classical solutions of the mPNP system was established.

Other forms of transformations though with inherent limitations to improve accuracy in the approximation of ion flow parameters in cell membranes have been suggested in several models such as steady state modified Poisson-Boltzmann (mPB) model. The mPB was later improved using Lambert-W special functions and the existence and uniqueness of its weak solution established. [5] equally derived a simple and effective modified Poisson Nernst Planck (mPNP) based on spartial modification of the diffusion coefficients of ions.

The work presented here considers a mathematical model by Lennard Jones (LJ) in 2D. The model consists of nonlinear PDEs with transport properties assumed to vary in continuous and differentiable manner. Critically in this study we employ the LJ repulsive potential for the finite size effects to device a further improved PNP equation. We use the variational approach to derive the total energy for LJ repulsive potential which leads to generation of a system of equations that incorporates contribution of finite size effects. Consequently, analysis of the local existence of weak solutions of the resultant mPNP by constructing an approximate solution in a finite dimensional space in L² is carried out.

2. Model Description

Deterministic mathematical model for simulating ion transport in bio cells are space and time dependent nonlinear PDEs posed as Initial value problems(IVP). Upon incorporating realistic and necessary physics for the flow phenomena the mPNP-IVP becomes complex deterring analytical determination of solution. In principle, to reduce intricacy in the approximation when fundamentally retaining accuracy in the estimation of the flow parameters, interaction of charged particles alone is declared sufficient. This leads to dropping inclusion of charges interaction with the channel walls and fluid which are assumed to have minimal contribution in the approximation results.

2.1. Modified Poisson Nernst-Planck (mPNP) Equation

The integral form of energy equation that integrates finite size effects and interactions between charged particles can be modelled as repulsive or attractive spherical particles in the energy term. The energy of these effects in the microscopic scale are summed to represent potential between the positive and negative charges [6] [7] [8] given by

$\begin{array}{c}E=\frac{1}{2}{\displaystyle \int }{\displaystyle {\int }_{\Omega }}{\Psi }_{nn}\left(\left|x-y\right|\right)c_n\left(x\right)c_n\left(y\right)\text{d}x\text{d}y\\ +\frac{1}{2}{\displaystyle \int }{\displaystyle {\int }_{\Omega }}{\Psi }_{np}\left(\left|x-y\right|\right)c_n\left(x\right)c_p\left(y\right)\text{d}x\text{d}y\\ +\frac{1}{2}{\displaystyle \int }{\displaystyle {\int }_{\Omega }}{\Psi }_{pp}\left(\left|x-y\right|\right)c_p\left(x\right)c_p\left(y\right)\text{d}x\text{d}y\end{array}$ (1)

where the repulsion or attractive potential between two balls of ions i and j of radius $a_i$, $a_j$ respectively situated at $x$ and $y$ in the two-dimensional spatial domain $\Omega$ given by

${\Psi }_{ij}\left(\left|x-y\right|\right)={\epsilon }_{ij}\frac{{\left(a_i+a_j\right)}^{12}}{{\left|x-y\right|}^{12}}\mathrm{,}$ (2)

${\epsilon }_{ij}$ is the chosen energy coupling constant and $c_n$, $c_p$ are the negative and positive charge densities, respectively.

Variation of Equation (1) is achieved by differentiating it with respect to charge density to result into flux due to the finite size effects of the charge densities, [7]. When the fluxes are added to the Nernst Planck (NP) equation results into modified time dependent equations that incorporates finite size effects each representing rates of time change in concentration of the negative and positive ions, respectively as in Equations (3) and (4) below.

$\begin{array}{c}\frac{\partial c_n}{\partial t}=\nabla \cdot \left[D_n\left(\nabla c_n-\frac{e}{K_{B}T}{c}_n\nabla \Phi \right)\right]\\ +\nabla \cdot [\frac{c_n\left(x\right)}{K_{B}T}{\displaystyle {\int }_{\Omega }}12{\epsilon }_{nn}\frac{{\left(a_n+a_n\right)}^{12}}{{\left|x-y\right|}^{14}}\left(x-y\right)c_n\left(y\right)\text{d}y\\ +{\displaystyle {\int }_{\Omega }}6{\epsilon }_{np}\frac{{\left(a_n+a_p\right)}^{12}}{{\left|x-y\right|}^{14}}\left(x-y\right)c_p\left(y\right)\text{d}y]\end{array}$ (3)

$\begin{array}{c}\frac{\partial c_p}{\partial t}=\nabla \cdot \left[D_p\left(\nabla c_p+\frac{e}{K_{B}T}{c}_p\nabla \Phi \right)\right]\\ +\nabla \cdot [\frac{c_p\left(x\right)}{K_{B}T}{\displaystyle {\int }_{\Omega }}12{\epsilon }_{pp}\frac{{\left(a_p+a_p\right)}^{12}}{{\left|x-y\right|}^{14}}\left(x-y\right)c_p\left(y\right)\text{d}y\\ +{\displaystyle {\int }_{\Omega }}6{\epsilon }_{np}\frac{{\left(a_n+a_p\right)}^{12}}{{\left|x-y\right|}^{14}}\left(x-y\right)c_p\left(y\right)\text{d}y]\end{array}$ (4)

Coupling Equations (3) and (4) with Poisson Equation (5) below we obtain the governing equation

$\nabla \cdot \left(ϵ\nabla \Phi \right)={\rho }_0+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{z}_{i}e{c}_i,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\forall i=n,p,$ (5)

in which $K_{B}T$ is the thermal energy, with $K_B$ the Boltzmann constant and T is the absolute temperature, ${\rho }_0$ is the charge density of the protein, $c_{p\mathrm{,}n}$ and $D_{p\mathrm{,}n}$ are the charge densities and diffusion coefficients for p and n ions respectively, $\Phi$ is the electrostatic potential, $ϵ$ is the dielectric coefficient, e is the unit charge, $z_i$ is the valence and N is the number of ions. The energy variational approach applied to the energy of LJ repulsive spheres ensures that the resulting mPNP system satisfies the energy dissipation law given by;

$\begin{array}{l}\frac{\partial }{\partial t}{\displaystyle \int }\left[K_{B}T{\displaystyle \underset{i=n,p}{\sum }}{c}_i\log c_i+\frac{1}{2}\left({\rho }_0+{\displaystyle \underset{i=n,p}{\sum }}{z}_{i}e{c}_i\right)\Phi +{\displaystyle \underset{i=n,p}{\sum }}\frac{c_i}{2}{\displaystyle \int }{\Phi }_{ij}{c}_j\text{d}y\right]\text{d}x\\ =-{\displaystyle \int }{\displaystyle \underset{i=n,p}{\sum }}\frac{D_i{c}_i}{K_{B}T}{\left|K_{B}T\frac{\nabla c_i}{c_i}+z_{i}e\nabla \Phi -{\displaystyle \underset{ij=n,p}{\sum }}\nabla {\displaystyle \int }{\Phi }_{ij}{c}_j\text{d}y\right|}^2\text{d}x\end{array}$ (6)

For brevity and without loss of generality, we assume non-standard values of the parameters to be $ϵ=e=K_B=T=D=1$, $z_n=-1$ and $z_p=1$, ${\rho }_0=0$, and the flow to be two dimensional with a unit thickness. This simplification enables linearization of Equations (3), (4) and (5) to obtain governing equation of the flow through the convergent-divergent ion channel given by

$\frac{\partial c_n}{\partial t}-\nabla \cdot \left[\nabla c_n-c_n\nabla \Phi \right]=-\nabla \cdot \left(c_{n}w\right)$ (7a)

$\frac{\partial c_p}{\partial t}-\nabla \cdot \left[\nabla c_p+c_p\nabla \Phi \right]=-\nabla \cdot \left(c_{p}v\right)\Omega \times \left(0,T\right),$ (7b)

$\Delta \Phi =\left(c_n-c_p\right)$ (7c)

where

$\begin{array}{c}v={\displaystyle {\int }_{\Omega }}12{\epsilon }_{n,n}\frac{{\left(a_n+a_n\right)}^{12}}{{\left|\left|x-y\right|\right|}^{14}}\left(x-y\right)c_n\left(y\right)\text{d}y\\ +{\displaystyle {\int }_{\Omega }}6{\epsilon }_{n,p}\frac{{\left(a_n+a_p\right)}^{12}}{{\left|\left|x-y\right|\right|}^{14}}\left(x-y\right)c_p\left(y\right)\text{d}y,\end{array}$

$\begin{array}{c}w={\displaystyle {\int }_{\Omega }}12{\epsilon }_{p,p}\frac{{\left(a_p+a_p\right)}^{12}}{{\left|\left|x-y\right|\right|}^{14}}\left(x-y\right)c_p\left(y\right)\text{d}y\\ +{\displaystyle {\int }_{\Omega }}6{\epsilon }_{n,p}\frac{{\left(a_n+a_p\right)}^{12}}{{\left|\left|x-y\right|\right|}^{14}}\left(x-y\right)c_n\left(y\right)\text{d}y.\end{array}$

2.2. Boundary and Initial Conditions

Consider a bounded domain $\Omega \subset R^d$, for $d\in \aleph$, we seek the solutions $\left(c_n\mathrm{,}{c}_p\mathrm{,}\Phi \right)\left(x\mathrm{,}t\right)$ in $\Omega$ to the Cauchy problem described by Equation (7) subject to initial data

$\left(c_n,\mathrm{<mi>\; </mi>}{c}_p,\mathrm{<mi>\; </mi>}\Phi \right)\left(x,0\right)=\left(c_{n0},\mathrm{<mi>\; </mi>}{c}_{p0},\mathrm{<mi>\; </mi>}{\Phi }_0\right)\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{for}x\in R^2$ (8)

We specify the Dirichlet boundary conditions to represent fixed electrostatic potential at the boundaries as;

$\Phi \left(0,t\right)=0,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\Phi \left(1,t\right)=1$ (9)

and prescribe Neumann boundary condition describing null charge density fluxes at the boundaries by;

$\frac{\Phi }{\partial \nu }=\frac{\partial c_i}{\partial \nu }=0$ (10)

where $i=n,p$; and $\nu$ is the unit outward normal. Throughout the paper we assume electro-neutrality conditions at the boundaries, implying that the;

$\underset{x\to \infty }{\lim }{c}_n\left(x,0\right)=c_p\left(x,0\right)$ (11)

3. Existence of Approximate Solution of mPNP

In this section, energy method is used to prove existence of solution of the governing equation for ion transport through cell membrane. This will start by first defining the space in which the solution is estimated, describing the local existence, determining prior energy estimate and finally working out the discrete solution and its convergence to the defined bound.

3.1. Local Existence of Solution

The approach involves determining the priori estimates on Sobolev norms of concentration $c_n$, $c_p$. Galerkin method, see [9] [10] is introduced to approximate the solution of mPNP equation by projection of the equation into finite dimensional subspace, $b_k$. The fundamental objective of this method in the proof of existence is to approximate $c_n\mathrm{,}{c}_p\in L^2\left(\mathrm{0,}T\right)$ by functions $c_{nk}\mathrm{,}{c}_{pk}\in b_k\left(\mathrm{0,}T\right)$ which takes values in a finite dimensional subspace $b_k\in L^2\left(\Omega \right)$ of dimension k.

We project Equation (7) onto $b_k$ to obtain $c_{nk}$, $c_{pk}$ which satisfies the equation upto a residual orthogonal to $b_k$. This gives rise to a system of Ordinary Differential Equations (ODEs) in $c_{nk}$, $c_{pk}$ which has a solution by standard ODE theory. The resultant solution $c_{nk}$ and $c_{pk}$ satisfies an energy estimate of the same form as a prior estimate for the solution of Equation (7). These estimates are uniform in $b_k$, and permits us to impose the limit $k\to \infty$ to obtain solution of Equation (7) in a bounded domain $\Omega \subset R^N$ for $N\le 2$. Sobolev spaces $H^k=w^{k\mathrm{,2}}\left(\Omega \right)$, $c\left(\Omega \right)$ denotes spaces of continuously differentiable and $c^k\left(\Omega \right)$ denotes k times continuously differentiable functions,in addition ${\Vert \mathrm{.}\Vert }_2$ denotes $L^2$ norm, ${\Vert \mathrm{.}\Vert }_1$ denotes $H^1$ norm and ${\Vert \mathrm{.}\Vert }_{\infty }$ denotes $L^{\infty }$ norm.

3.2. The Priori Energy Estimate

In general, it is demanding to solve mPNP analytically because of nonlinearity, thus derivation of some energy estimates for the solutions of the system of Equation (7) by assuming v and w are given functions becomes a possible way for studying and analysing the physical problem.

Lemma 1 According to [3], given $\Omega \subset R^2$, taking $T>0$ and $v\mathrm{,}w\in L^2\left(\mathrm{0,}T\mathrm{<mo>\; ;\; </mo>}{L}^2\right)$, then if $\left(c_n\mathrm{,}{c}_p\right)\in L^{\infty }\left(\mathrm{0,}T\mathrm{<mo>\; ;\; </mo>}{L}^2\right)\cap L^2\left(\mathrm{0,}T\mathrm{<mo>\; ;\; </mo>}{H}_0^1\right)$ there exists a constant $\beta >0$, such that for any $t\in \left[\mathrm{0,}T\right]$

$\begin{array}{l}{\Vert \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }_2^2+{\displaystyle {\int }_0^t}{\Vert \nabla \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \\ \le ex{p}^{2{\Vert \nabla \Phi \Vert }_2^{2}t}\left[{\Vert \left(c_{n0}\mathrm{,}{c}_{p0}\right)\left(x\mathrm{,0}\right)\Vert }_2^2+{\displaystyle {\int }_0^t}\left(\beta {\Vert \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }_1+D{\Vert \left(w\mathrm{,}v\right)\left(x\mathrm{,}\tau \right)\Vert }_1\text{d}\tau \right)\right]\mathrm{,}\end{array}$

and

$\begin{array}{l}{\Vert \nabla \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }_2^2+\gamma {\displaystyle {\int }_0^t}{\Vert \Delta \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \\ \le ex{p}^{Kt}\left[{\Vert \nabla \left(c_{n0}\mathrm{,}{c}_{po}\right)\Vert }_2^2+D{\displaystyle {\int }_0^t}{\Vert \nabla \left(w\mathrm{,}v\right)\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \right]\mathrm{.}\end{array}$

Proof Multiplying Equation (7a) by $c_n$ and integrate over the domain $\Omega$, then using young’s inequality to simplify we obtain

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert c_n\Vert }_2^2+{\Vert \nabla c_n\Vert }_2^2+{\displaystyle {\int }_{\Omega }}\nabla \cdot \left(c_n\nabla \Phi \right)c_n=-{\displaystyle {\int }_{\Omega }}\nabla \cdot \left(c_{n}w\right)c_n\mathrm{.}$ (12)

Integrating by parts and applying holder’s inequality, ${\displaystyle {\int }_{\Omega }}\left|ab\right|\le \Vert a\Vert \Vert b\Vert$, in the last two terms of Equation (12) we get

$(\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left|\nabla \cdot \left(c_n\nabla \Phi \right)c_n\right|\le {\Vert c_n\Vert }_2{\Vert \nabla \Phi \Vert }_2{\Vert \nabla c_n\Vert }_2\\ -{\displaystyle {\int }_{\Omega }}\left|\nabla \cdot \left(c_{n}w\right)c_n\right|\le {\Vert c_n\Vert }_2{\Vert w\Vert }_2{\Vert \nabla c_n\Vert }_2\end{array}$ (13)

then by substituting in Equation (13) followed by imposition of Cauchy’s inequality in Equation (12) results into

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert c_n\Vert }_2^2+{\Vert \nabla c_n\Vert }_2^2\le \frac{1}{4\xi }{\Vert c_n\Vert }_2^2{\Vert \nabla \Phi \Vert }_2^2+\xi {\Vert \nabla c_n\Vert }_2^2+\frac{1}{4\xi }{\Vert c_{n}w\Vert }_2^2+\xi {\Vert \nabla c_n\Vert }_2^2$ (14)

by further substituting $\xi =\frac{1}{4}$ to simplify right hand side of Equation (14) we finally get

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert c_n\Vert }_2^2+\frac{1}{2}{\Vert \nabla c_n\Vert }_2^2\le {\Vert c_n\Vert }_2^2{\Vert \nabla \Phi \Vert }_2^2+\xi {\Vert c_n\Vert }^4+D_{\xi }{\Vert w\Vert }^4\mathrm{.}$ (15)

Similarly multiplying Equation (7b) by $c_p$ and integrating by parts over the domain we obtain equivalence of Equation (12) as below

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert c_p\Vert }_2^2+{\Vert \nabla c_p\Vert }_2^2-{\displaystyle {\int }_{\Omega }}\nabla \cdot \left(c_p\nabla \Phi \right)c_p=-{\displaystyle {\int }_{\Omega }}\nabla \cdot \left(c_{p}v\right)c_p$ (16)

then following the same procedure as Equations (13) and (14) in Equation (16) we obtain

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert c_p\Vert }_2^2+\frac{1}{2}{\Vert \nabla c_p\Vert }_2^2\le {\Vert c_p\Vert }_2^2{\Vert \nabla \Phi \Vert }_2^2+\xi {\Vert c_p\Vert }^4+D_{\xi }{\Vert v\Vert }^4\mathrm{.}$ (17)

Applying Gronwall inequality in Equations (15) and (17) when taking $\beta >0$ and $C>0$ gives

$\begin{array}{l}{\Vert \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }_2^2+{\displaystyle {\int }_0^t}{\Vert \nabla \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \\ \le ex{p}^{2{\Vert \nabla \Phi \Vert }_2^{2}t}\left[{\Vert \left(c_{n0}\mathrm{,}{c}_{p0}\right)\left(x\mathrm{,0}\right)\Vert }_2^2+{\displaystyle {\int }_0^t}\left(\beta {\Vert \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }^4+D{\displaystyle {\int }_0^t}{\Vert \left(w\mathrm{,}v\right)\left(x\mathrm{,}\tau \right)\Vert }^4\text{d}\tau \right)\right]\end{array}$ (18)

then by continuous embedding of $H_1$ into $L_4$, Equation (18) is simplified to obtain

$\begin{array}{l}{\Vert \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }_2^2+{\displaystyle {\int }_0^t}{\Vert \nabla \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \\ \le ex{p}^{2{\Vert \nabla \Phi \Vert }_2^{2}t}\left[{\Vert \left(c_{n0}\mathrm{,}{c}_{p0}\right)\left(x\mathrm{,0}\right)\Vert }_2^2+{\displaystyle {\int }_0^t}\left(\beta {\Vert \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }_1+D{\Vert \left(w\mathrm{,}v\right)\left(x\mathrm{,}\tau \right)\Vert }_1\text{d}\tau \right)\right]\mathrm{.}\end{array}$ (19)

To obtain the energy estimates for the derivatives we multiply Equation (7a) by $\Delta c_n$ and integrate by parts over the spatial domain $\Omega$ we obtain

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert \nabla c_n\Vert }_2^2+{\Vert \Delta c_n\Vert }_2^2+{\displaystyle {\int }_{\Omega }}\nabla \cdot \left(c_n\nabla \Phi \right)\Delta c_n=-{\displaystyle {\int }_{\Omega }}\nabla \cdot \left(c_{n}w\right)\Delta c_n\mathrm{.}$ (20)

Re-evaluating two last terms of Equation (20) using Cauchy’s inequality give rise to Equations (21) and (22) below:

$\begin{array}{c}{\displaystyle \int }\nabla \cdot \left(c_n\nabla \Phi \right)\Delta c_n=\Delta c_n\left[c_n\Delta \Phi +\left(\nabla \cdot c_n\right)\nabla \Phi \right]\\ \le {\Vert \Delta c_n\Vert }_2\left[{\Vert c_n\Vert }_{\infty }{\Vert \Delta \Phi \Vert }_2+{\Vert \nabla \Phi \Vert }_{\infty }{\Vert \left(\nabla \cdot c_n\right)\Vert }_2\right]\\ \le D_1{\Vert \Delta c_n\Vert }_2{\Vert \Delta \Phi \Vert }_2+D_2{\Vert \Delta c_n\Vert }_2\Vert \nabla c_n\Vert \\ \le D_1{\Vert \Delta c_n\Vert }_2{\Vert \Delta \Phi \Vert }_2+D_2\left[\xi {\Vert \nabla c_n\Vert }_2^2+\frac{1}{4\xi }{\Vert \Delta c_n\Vert }_2^2\right]\end{array}$ (21)

$\begin{array}{c}{\displaystyle \int }\nabla \cdot \left(c_{n}w\right)\Delta c_n=\Delta c_n\left[c_n\nabla w+\left(\nabla \cdot c_n\right)w\right]\\ \le {\Vert \Delta c_n\Vert }_2\left[{\Vert c_n\Vert }_{\infty }{\Vert \nabla w\Vert }_2+{\Vert w\Vert }_{\infty }{\Vert \left(\nabla \cdot c_n\right)\Vert }_2\right]\\ \le D_3{\Vert \Delta c_n\Vert }_2{\Vert \nabla w\Vert }_2+D_4{\Vert \Delta c_n\Vert }_2{\Vert \nabla c_n\Vert }_2\\ \le D_3\left[\xi {\Vert \Delta c_n\Vert }_2^2+\frac{1}{4\xi }{\Vert \nabla w\Vert }_2^2\right]+D_4\left[\xi {\Vert \nabla c_n\Vert }_2^2+\frac{1}{4\xi }{\Vert \Delta c_n\Vert }_2^2\right]\end{array}$ (22)

substituting Equations (21) and (22) in Equation (20) we have

$\begin{array}{l}\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert \nabla c_n\Vert }_2^2+{\Vert \Delta c_n\Vert }_2^2\\ \le D_1{\Vert \Delta c_n\Vert }_2{\Vert \Delta \Phi \Vert }_2+D_2\left[\xi {\Vert \nabla c_n\Vert }_2^2+\frac{1}{4\xi }{\Vert \Delta c_n\Vert }_2^2\right]\\ +D_3\left[\xi {\Vert \Delta c_n\Vert }_2^2+\frac{1}{4\xi }{\Vert \nabla w\Vert }_2^2\right]+D_4\left[\xi {\Vert \nabla c_n\Vert }_2^2+\frac{1}{4\xi }{\Vert \Delta c_n\Vert }_2^2\right]\mathrm{,}\end{array}$ (23)

and upon further simplification Equation (23) give rise to (24)

$\frac{\text{d}}{\text{d}t}{\Vert \nabla c_n\Vert }_2^2+k_1\tilde{\xi }{\Vert \Delta c_n\Vert }_2^2\le k_2\xi {\Vert \nabla c_n\Vert }_2^2+k_3{\Vert \nabla w\Vert }_2^2+k_4{\Vert \Delta c_n\Vert }_2{\Vert \Delta \Phi \Vert }_2\mathrm{.}$ (24)

Similarly multiplying Equation (7b) by $\Delta c_p$ and following the same procedure as in Equations (21)-(23) we get

$\frac{\text{d}}{\text{d}t}{\Vert \nabla c_p\Vert }_2^2+k_5\tilde{\xi }{\Vert \Delta c_p\Vert }_2^2\le k_6\xi {\Vert \nabla c_p\Vert }_2^2+k_7{\Vert \nabla v\Vert }_2^2+k_8{\Vert \Delta c_p\Vert }_2{\Vert \Delta \Phi \Vert }_2$ (25)

then adding Equations (24) and (25) we have

$\frac{\text{d}}{\text{d}t}{\Vert \nabla \left(c_n\mathrm{,}{c}_p\right)\Vert }_2^2+\gamma {\Vert \Delta \left(c_n\mathrm{,}{c}_p\right)\Vert }_2^2\le K{\Vert \nabla \left(c_n\mathrm{,}{c}_p\right)\Vert }_2^2+D{\Vert \nabla \left(w\mathrm{,}v\right)\Vert }_2^2$ (26)

where the constants $\gamma ,D>0$ and $K>0$.

Finally when Gronwall inequality is imposed on Equation (26) we obtain

$\begin{array}{l}{\Vert \nabla \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}t\right)\Vert }_2^2+\gamma {\displaystyle {\int }_0^t}{\Vert \Delta \left(c_n\mathrm{,}{c}_p\right)\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \\ \le ex{p}^{Kt}\left[{\Vert \nabla \left(c_{n0}\mathrm{,}{c}_{p0}\right)\Vert }_2^2+D{\displaystyle {\int }_0^t}{\Vert \nabla \left(w\mathrm{,}v\right)\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \right]\mathrm{.}\end{array}$ (27)

ending proof of the lemma.

3.3. Galerkin Approximation

Next we use an ODE theory to get local in time solution to a finite dimensional approximation $c_{nk}$, $c_{pk}$ of Equation (7) satisfying the same energy bounds as $c_n$, $c_p$. Using Galerkin approximation we first construct finite dimensional approximation to Equation (7). Let $e_j=e_j\left(x\right)$ be smooth function and ${\left\{e_j\right\}}_{j=1}^{\infty }$ be an orthonormal basis for the functions in $L^2\left(\Omega \right)$. Let $b_k\in L^2\left(\Omega \right)$ be a finite dimensional subspace given by

$b_k=\left\{u\in L^2\left(\Omega \right);u\left(x\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}}{d}_j{e}_j\left(x\right)\right\}.$ (28)

We define a projection operator $q_n$ from $L^2\left(\Omega \right)$ to $b_k$ and if $u={\displaystyle {\sum }_{j=1}^{\infty }}{u}_j\left(t\right)e_j\left(x\right)$, then $q_{n}u={\displaystyle {\sum }_{j=1}^n}{u}_j\left(t\right)e_j\left(x\right)$ for $u\in L^2\left(\Omega \right)$. From the properties of orthogonal projection operator $\Vert q_{n}u\Vert \le \Vert u\Vert$ holds for $u\in b_k$.

Lemma 2 For any $m>0$ and $z\in L^2\left(\Omega \right)$,

$\left({\nabla }_{m}u\mathrm{,}{\nabla }_m{q}_{n}z\right)=\left({\nabla }_{m}u\mathrm{,}{\nabla }_{m}z\right)$, $\forall u\in b_k$.

Lemma 3 If $z_n\to z$ in $c^1\left(\Omega \right)$, then ${\Vert q_n{z}_n-z\Vert }_{\infty }\to 0$ as $n\to \infty$.

The proof of Lemma 2 and 3 is trivial see, [11].

Taking $s=c_n,c_p$; $s_n=c_{nk},c_{pk}$; and $p=\left(v\mathrm{,}w\right)$ in Equations (7a) and (7b), we define Galerkin approximation by seeking solution $s_n\in b_k$ with finite representation given by

$s_n\left(x,t\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}}{a}_j\left(t\right)e_j\left(x\right)$ (29)

satisfying

${\partial }_t{s}_n-\Delta s_n=\nabla \cdot q_n\left[s_n\nabla \Phi \right]-\nabla \cdot q_n\left[s_{n}p\right]$ (30)

subject to initial conditions

$s_n\left(x,0\right)=q_n{s}_0.$ (31)

Through transformation we obtain equivalence of Equation (30) given by

$\frac{\text{d}}{\text{d}t}{a}_j-{e^{″}}_j{a}_j={e^{\prime }}_j\left({\Phi }^{\prime }-w\right)+e_j{a}_j\left[{\Phi }^{″}-w^{\prime }\right]$ (32)

with initial conditions given by

$a_j\left(0\right)=\left(e_j{s}_0\right).$ (33)

The integral form of Equations (30) and (31) becomes

$s_n\left(x\mathrm{,}t\right)=q_n{s}_0+{\displaystyle {\int }_0^t}\left\{\Delta s_n+\nabla \cdot q_n\left[s_n\nabla \Phi -s_{n}p\right]\left(x\mathrm{,}\tau \right)\right\}\text{d}\tau$ (34)

Theorem 1 For any integer n and real number $T>0$ there exists a unique solution $s_n$ to partial differential Equations (30) subject to initial condition (31) for $t\in \left(\mathrm{0,}T\right)$ also satisfying the ordinary differential Equations (32) with initial condition given by Equation (33).

Proof Considering equivalence of Equations (30) and (32) subject to respective initial conditions as above, proof of existence of local solution when taking any fixed $n\ge 1$ for the system of Equations (30)-(31) in the interval $t\in \left(\mathrm{0,}T\right)$ to imply existence of solution to the Equations (32)-(33). Hence, we proceed to multiply Equation (32) by $s_n$ and integrate by parts in the domain to obtain

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert s_n\Vert }^2+{\Vert \nabla s_n\Vert }^2={\displaystyle {\int }_{\Omega }}\nabla \cdot q_n\left[s_n\nabla \Phi \right]\nabla s_n+{\displaystyle {\int }_{\Omega }}\nabla \cdot q_n\left[s_{n}p\right]\nabla s_n\mathrm{.}$ (35)

Following integration process as in Equation (13), the last two terms of Equation (35) reduces to

$\{\begin{array}{l}{\displaystyle {\int }_{\Omega }}{q}_n\left[s_n\nabla \Phi \right]\nabla s_n={\displaystyle {\int }_{\Omega }}\left[s_n\nabla \Phi \right]\nabla s_n={\Vert s_n\Vert }_2{\Vert \nabla \Phi \Vert }_2{\Vert \nabla s_n\Vert }_2,\text{and}\\ {\displaystyle {\int }_{\Omega }}{q}_n\left[s_{n}p\right]\nabla s_n={\displaystyle {\int }_{\Omega }}\left[s_{n}p\right]\nabla s_n={\Vert s_n\Vert }_2{\Vert p\Vert }_2{\Vert \nabla s_n\Vert }_2.\end{array}$ (36)

Substituting Equation (36) in (35) then using Cauchy’s inequality gives

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert s_n\Vert }^2+{\Vert \nabla s_n\Vert }^2\le \frac{1}{4\xi }{\Vert s_n\Vert }^2{\Vert \nabla \Phi \Vert }^2+\xi {\Vert \nabla s_n\Vert }^2+\frac{1}{4\xi }{\Vert s_n\Vert }^2{\Vert w\Vert }^2+\xi {\Vert \nabla s_n\Vert }^2$ (37)

then taking $\xi =\frac{1}{4}$ and applying continuous embedding of $H_1$ into $L_4$, [9] Equation (36) becomes

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{\Vert s_n\Vert }^2+{\Vert \nabla s_n\Vert }^2\le {\Vert s_n\Vert }^2{\Vert \nabla \Phi \Vert }^2+{\Vert s_n\Vert }^2{\Vert p\Vert }^2\\ \le {\Vert s_n\Vert }^2{\Vert \nabla \Phi \Vert }^2+{\Vert s_n\Vert }^4+{\Vert p\Vert }^4\\ \le {\Vert s_n\Vert }^2{\Vert \nabla \Phi \Vert }^2+{\Vert s_n\Vert }_1+{\Vert p\Vert }_1\mathrm{.}\end{array}$ (38)

Consequently applying Gronwall inequality in Equation (38) gives

${\Vert s_n\Vert }^2+{\displaystyle {\int }_0^t}{\Vert \nabla s_n\Vert }^2\text{d}\tau \le ex{p}^{2{\Vert \nabla \Phi \Vert }^{2}t}\left[{\Vert s_0\Vert }^2+{\displaystyle {\int }_0^t}{\Vert s_n\left(\tau \right)\Vert }_1+{\Vert p\Vert }_1\left(\tau \right)\text{d}\tau \right]$ (39)

and from parseval’s inequality given in [9],

${\displaystyle \underset{j=1}{\overset{n}{\sum }}}{a}_j^2\le \infty$ (40)

we finally have $a_j\left(t\right)$ that can be determined for any $t\in \left(\mathrm{0,}T\right)$.

Theorem 2 For an integer n, the solution $s_n$ in Theorem 1 satisfies the same priori energy bounds as on s.

Proof; Using Lemma 2 to eliminate $q_n$ in the inner product of Equation (30) and similar procedure as in Equations (20)-(27) to find the time integral estimate of $s_n$, we obtain the bounds of derivatives given by

${\Vert \nabla s_n\left(x\mathrm{,}t\right)\Vert }_2^2+\hat{\gamma }{\displaystyle {\int }_0^t}{\Vert \Delta s_n\left(x\mathrm{,}\tau \right)\Vert }_2^2\text{d}\tau \le ex{p}^{\alpha t}\left[{\Vert \nabla s_0\Vert }_2^2+C{\displaystyle {\int }_0^t}{\Vert \nabla p\Vert }_2^2\text{d}\tau \right]\mathrm{.}$ (41)

To prove boundedness of the derivatives in $s_n$ we introduce the following lemma;

Lemma 4 For an integer $n>1$ and $t\in \left(\mathrm{0,}T\right)$, the solution $s_n$ in Theorem 1 also satisfies

$\Vert {\partial }_t{s}_n\Vert \mathrm{<mo>\; ;\; </mo>}\Vert {\partial }_t\nabla s_n\Vert \le E$ (42)

with E independent of n and T.

Proof Consider $L^2$ estimate of Equation (30) and noting that $q_n$ vanishes with respect to Lemma 2, we get

$\begin{array}{c}{\Vert {\partial }_t{s}_n\Vert }_2\le {\Vert \Delta s_n\Vert }_2+{\Vert \nabla \cdot \left(s_n\nabla \Phi \right)\Vert }_2+{\Vert \nabla \cdot \left(s_{n}p\right)\Vert }_2\\ \le {\Vert \Delta s_n\Vert }_2+{\Vert \nabla s_n\Vert }_2{\Vert \nabla \Phi \Vert }_{\infty }+{\Vert s_n\Vert }_{\infty }{\Vert \Delta \Phi \Vert }_2+{\Vert \nabla s_n\Vert }_2{\Vert p\Vert }_{\infty }+{\Vert s_n\Vert }_{\infty }{\Vert \nabla p\Vert }_2\\ \le E\end{array}$ (43)

Evaluating the derivative of Equation (43) with respect to x gives

$\begin{array}{c}{\Vert {\partial }_t\nabla s_n\Vert }_2\le {\Vert {\nabla }^3{s}_n\Vert }_2+2{\Vert \nabla s_n\Vert }_2{\Vert \Delta \Phi \Vert }_{\infty }+{\Vert s_n\Vert }_{\infty }{\Vert {\nabla }^3\Phi \Vert }_2+{\Vert \Delta s_n\Vert }_2{\Vert \nabla \Phi \Vert }_{\infty }\\ +2{\Vert \nabla s_n\Vert }_2{\Vert \nabla p\Vert }_{\infty }+{\Vert s_n\Vert }_{\infty }{\Vert \Delta p\Vert }_2+{\Vert \Delta s_n\Vert }_2{\Vert p\Vert }_{\infty }\\ \le E\end{array}$ (44)

Therefore the derivatives of $s_n$ are uniformly bounded

It is imperative to remark that a uniform bounds on $s_n$ leads to uniform bounds on the derivatives. This uniform boundedness proves equicontinuity of $s_n$. Thus by Arzela Ascoli theorem [12] there exists a subsequence $s_{nj}$ which converges uniformly to s on the interval $t\in \left(\mathrm{0,}{T}_n\right)$.

Theorem 3 Let $s_n$ is a sequence of functions $R^n\to R$ which is equicontinous and bounded. Then there exists a subsequence ${\left\{s_{nj}\right\}}_{j=1}^{\infty }$ which converges in the sup norm to s, where s satisfies the integral equation

$s\left(x\mathrm{,}t\right)=s_0\left(x\right)+{\displaystyle {\int }_0^t}\left\{\Delta s\left(x\mathrm{,}\tau \right)-\nabla \cdot \left(s\nabla \Phi -sp\right)\left(x\mathrm{,}\tau \right)\right\}\text{d}\tau$ (45)

implying that s also satisfies Equation (7).

Proof Since $s_n$ and its derivatives have bounds independent of t and n, it follows that ${\left\{s_n\right\}}_{n=1}^{\infty }$ forms an equicontinous family of functions which have a

subsequence ${\left\{s_{nj}\right\}}_{j=1}^{\infty }$ that converges in the sup norm for each $t\in \left(\mathrm{0,}T\right)$ to s. Since $s_{nj}\to s\in C\left(\left(\mathrm{0,}T\right)\mathrm{,}{L}^2\right)$ it follows that ${\Vert s_{nj}\nabla \Phi -s_n\nabla \Phi \Vert }_{\infty }\to 0$ as $j\to \infty$, ${\Vert \Delta s_{nj}-\Delta s_n\Vert }_{\infty }$ as $j\to \infty$. Therefore each time integral

${\displaystyle {\int }_0^t}\left\{\Delta s_{nj}+\nabla \cdot q_{nj}\left[s_{nj}\nabla \Phi -s_{nj}p\right]\left(x\mathrm{,}\tau \right)\right\}\text{d}\tau$ (46)

is bounded independent in t. For any given x, ${\left\{s_{nj}\right\}}_{j=1}^{\infty }$ converges uniformly for $t\in \left[\mathrm{0,}T\right]$ as $j\to \infty$. Therefore from dominating convergence and using lemma 3, we obtain

$\begin{array}{l}\underset{j\to \infty }{lim}{\displaystyle {\int }_0^t}\left\{\Delta s_{nj}+\nabla \cdot q_{nj}\left[s_{nj}\nabla \Phi -s_{nj}p\right]\left(x\mathrm{,}\tau \right)\right\}\text{d}\tau \\ ={\displaystyle {\int }_0^t}\left\{\Delta s_n+\nabla \cdot \left[s_n\nabla \Phi -s_{n}p\right]\left(x\mathrm{,}t\right)\right\}\text{d}\tau \end{array}$ (47)

hence $li{m}_{j\to \infty }{s}_{nj}\left(x\mathrm{,}t\right)=s\left(x\mathrm{,}t\right)$ and s satisfies Equation (7).

Having proved the existence and convergence of approximate solution, its important to note that the uniqueness of this solution conforms to proposition 3.9 on uniqueness found in Equations (3.39)-(3.40), see [3].

4. Conclusion

In this study mPNP differential equations which inco-orporate volume size effects is derived. The system of differential equations linearized and the proof of existence of its approximate analytical solution done. Lastly, a projection operator is used to map differential equation from L² spaces into a finite dimensional subspace $b_k\in L^2\left(\Omega \right)$ of dimension k generating equations whose solutions converge to that of mPNP equations. In the next study we will examine the existence and uniqueness of numerical approximate solution using Galerkin approach for which we shall conduct a numerical experiment for a 2D flow through ion channels.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Karniadakis, G., Beskok, A. and Narayan, A. (2005) Microows and Nanoows: Fundamentals and Simulation. In: Antman, S.S., Marsden, J.E. and Sirovich, L., Eds., Springer-Verlag, Berlin, Vol. 29, 30-36.
[2]	Qiao, Y., Bin, T. and Benzhuo, L. (2014) Ion Size Effects to Molecular Solvation Energy and to Ion Current across a Channel Resulted from the Non-Uniform SIZE Modified PNP System. The Journal of Chemical Physics, 140, Article No. 174102. https://doi.org/10.1063/1.4872330
[3]	Yong, Y., Chun, L. and Zhang, T. (2016) A Generalized Poisson Nernst-Planck Navier Stokes Model on the UID with the Crowded Charged Particles; Derivation and Its Well Posedness. SIAM Journal on Mathematical Analysis, 48, 3191-3235. https://doi.org/10.1137/16M1055104
[4]	Hsieh, C., Yunkyong, H., Hijin, L., Tai-Chia, L. and Chun, L. (2014) Transport of Charged Particles: Entropy Production and Maximum Dissipation Principle. Journal of Mathematical Analysis and Applications, 422, 31.
[5]	Farjana, S., Zhongming, W. and Shenggao, Z. (2017) A Modified Poisson Nernst-Planck Model with Excluded Volume Effects; Theory and Numerical Implementation. Communications in Mathematical Sciences, 16, 26.
[6]	Eisenberg, B. (2011) Mass Action in Ionic Solutions. Chemical Physics Letters, 511, 1-6. https://doi.org/10.1016/j.cplett.2011.05.037
[7]	Eisenberg, B., Yunkyong, H. and Chun, L. (2010) Energy Variational Analysis of Ions in Water and Channels; Field Theory for Primitive Models of Complex Ionic UIDS. The Journal of Chemical Physics, 133, Article No. 104104. https://doi.org/10.1063/1.3476262
[8]	Yunkyong, H., Eisenberg, B. and Liu, C. (2011) A Mathematical Model for Hard Sphere Repulsion in Ionic Solutions. Communications in Mathematical Sciences, 9, 459-475. https://doi.org/10.4310/CMS.2011.v9.n2.a5
[9]	Hunter, J. (2014) Notes on Partial Differential Equations. 177-193.
[10]	Tomastz, D. (2008) Certain Methods of Proving Existence and Uniqueness in Partial Differential Equations Theory. Silesian University, Poland, 18.
[11]	Tanveer, S. (2017) Introduction to Energy Methods for Existence of Solution. Ohio state University, Columbus, 14.
[12]	Brezis, H. (2010) Functional Analysis, Sobolev Spaces and Partial Diffrential Equations. Springer, New York, NY, 201-312. https://doi.org/10.1007/978-0-387-70914-7_8