Solutions of Poisson-Nernst Planck Equations with Ion Interaction

Abidha Monica Gwecho¹, Wang Shu¹, Onyango Thomas Mboya², Sudheer Khan¹

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

DOI: 10.4236/am.2022.133020                

Abstract

This paper develops a mixed Finite Element Method (mFEM) based on both classical rectangular elements (with equal nodal points for all degrees of freedom) and Taylor-hood elements to solve Poisson-Nernst-Planck (PNP) equations with steric effects. The Nernst-Planck (NP) equation for ion fluxes is modified to incorporate finite-size effects in terms of hard-sphere repulsion. The resultant modified NP and Poisson equation representing electrostatic potential is then coupled to form a system of the equation that describes a realistic transport phenomenon in an ion channel. Consequently, we apply mFEM based on both Taylor-hood and classical rectangular elements to discretize the 2D steady system of equations with iterative linearization for the nonlinear components. The numerical scheme is first validated using a 2D analytical solution for PNP equations, the steric components varied and the effects on concentration analyzed and compared against PNP and modified PNP for two monovalent ion species. Classical rectangular elements presented a better comparable approximate solution than Taylor-hood.

Keywords

Modified PNP, Lennard Jonnes, Hard-Sphere, Taylor-Hood

1. Background

Biological ion channels are proteins located in the membranes of the cells of our body with a central pore that allows and controls the passage of charged ions such as calcium (Ca²⁺), potassium (K⁺), Sodium (Na²⁺) and so on. These channels support transfer across the membrane when opened by an appropriate stimulus. The concentration of these ions is one of the crucial factors in regulating many physiological functions including signal transfer in the nervous system, regulating the secretion of hormones in the cells, etc. Thus furthering understanding of the dynamics and mechanism of the ion channel is vital for medical intervention in a broad range of diseases areas and the determination of new research trends in medical fields [1].

One of the most popular continuum theories describing ion transport in such a complex biological system is the Poisson Nernst-Planck (PNP) model [2]. The PNP equations describe the electro diffusion of ions under the effect of electric field induced by ions charge themselves. The model uses a combination of, firstly, Poisson’s equation to describe how ions and the channel protein to create electrical potential and secondly, Nernst-Planck equations to describe migration and diffusion of ions in gradients of concentration and electrical potential. Furthermore, the two components of the PNP model form, also known as the “drift-diffusion equations” which is widely used to describe behaviour of semiconductors, solid-state devices like transistors, nano-devices in biophysics and physical chemistry as illustrated in [2] [3] [4] [5].

Among the popularly applied numerical techniques for PNP model include Finite Difference (FD) [5] [6] [7], spectral methods [4] [8] and Finite Element Methods (FEM) [9] [10] [11] [12] [13]. The most preferred and appropriate numerical method amongst them for biological channels is FEM, this is due to the inherent capability to adequately handle irregular geometries and non-uniform boundaries. A number of researchers have applied different types of FEM in the study of steric PNP models, for example, [9] generalized the Borukhov model to obtain a size modified PNP (SMPNP) that was able to treat non-uniform particle sizes. The objective was to show that: 1) Size Modified Poisson-Boltzmann (SMPB) model was implied by SMPNP equations under certain boundary and interface conditions and can be reproduced through the numerical solution of the SMPNP. 2) The side effects in the SMPNP effectively reduce the densities of highly concentrated counter-ions around the biomolecule. An accurate finite element method and convergent Gummel iteration were developed for the numerical solution of a completely coupled non-linear system of SMPNP.

The validity of the PNP model was extended by adding an Excess Chemical Potential (ECP) correction to account for finite ion size and water occupation by [4]. The modification of the standard drift-diffusion solution methodology was accomplished by adding an outer iteration for the correction, achieving feasible convergence rates with simple test structures. Under equilibrium conditions in the absence of fixed charge on the membrane, PNP theory predicted uniform ion densities while PNP/ECP predicted non-uniform charge distribution. [10] presented a mathematical model for finite-size effects using a regularized Lennard Jones (LJ) repulsive potential under an energy variational framework and its numerical verification to recover layering behaviour. Edge finite element method was used to solve a system of modified PNP and employed convex iteration scheme to ensure self-consistency between ionic concentrations and electrostatic potentials. Both PNP with LJ repulsive potential and Density function theory hard-sphere potential have the same overall behaviour of ion concentration. Ion concentration near the boundary was found to be larger than in bulk.

Thereafter, [13] proposed linearized local conservative mixed FEM for solving time-dependent PNP equations, where the NP and Poisson fluxes were introduced as new vector-valued variables. These methods were employed for spatial discretization, backward Euler with linearisation for temporal discretization hence the system decoupled and solved iteratively for each time step. Additionally, [14] developed a modified PNP macroscopic model accounting for size exclusion effects resulting in nonlinear mobilities. The numerical simulation of the model indicated decreasing behaviour in conductance which leads to current saturation, behaviour was not detected in the PNP case. Finally, it is important to note that FEM simulation of PNP with steric effects is mostly limited to 1D cases see [10]. This is despite realistic computationally expensive 3D ion channel see [6], a major contributing factor being difficulty in numerical convergence.

Consequently, this study adopts approximation of the LJ potential in the energy functional with hard-sphere repulsion using band limit function which cut-off higher frequencies and preserve spatial frequencies in Fourier modes. The approximation is done with help of a cutoff length (σ) approaching zero, see [8] to simplify the kernel in the functional. As a result, the improved energy function is used to derive equations of PNP with steric effects (mPNP), which is simple and numerically viable. The paper proposes and analyzes the mixed FEM method based on Taylor-hood elements and classical FEM using the bilinear element with equal degrees of freedom to discretize the domain. Furthermore, FEM algorithm is enforced over the domain with iterative linearization to solve the decoupled system of highly nonlinear mPNP and PNP equations for two monovalent ion species in 2D. Accordingly, the study analyses their concentration behaviour in the channel with changing steric values. In the previous study, we analyzed the existence of approximate solution of a 2D mPNP system, see [15] and numerical simulation using mFEM in [16], the current study intends to improve the accuracy of numerical approximate solutions obtained in the previous simulation and offer a validation scheme.

2. Mathematical Model

Given that classical PNP model have known limitations, for example, it neglects the finite size effects of ion particles and does not account for non-electrostatic interaction between ions see, [4] [9] [17] [18]. As a result the PNP model cannot adequately describe the flow behaviour of biological ion channel calling for improvement. Consequently, in recent years, mathematical studies in this area have modified the PNP equation to include finite volume effects and electrostatic interactions.

2.1. Poisson-Nernst-Planck Equation

The Poisson-Nernst-Planck mathematical model for analyzing ion transport is given by:

$\frac{\partial c_i}{\partial t}=\nabla \cdot D_i\left(\nabla c_i+\frac{z_{i}e}{K_{B}T}{c}_i\nabla \Phi \right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}i=1,\cdots ,N$ (2.1a)

$-\nabla \cdot \left(\epsilon \nabla \Phi \right)={\rho }_0+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{z}_{i}e{c}_i.$ (2.1b)

In the Nersnt-Planck (NP) Equation (2.1a), $c_i$ and $z_i$ are concentration and valence for the ith ion species. $\varphi$ is the electrostatic potential, $K_B$ is the Boltzmann constant, T is the absolute temperature, N is the number of ions, e is the unit charge and $D_i$ is the diffusion coefficient of ion species. In the Poisson Equation (2.1b), $\epsilon$ is the coupling energy constant and ${\rho }_0$ representing permanent charge density.

2.2. Poisson-Nernst-Planck Equation with Steric Effects

This subsection considers a continuum flow in a two dimensional problem to represent a channel in cell membrane. The energy variational approach is used to derive a system of differential equations including finite size effects of ions using the Lennard Jonnes repulsive potential. This potential introduces the ion interaction which is modeled as hard spheres in a rectangular space with unit thickness. The contribution of the repulsive potential to the total free energy functional is given by:

$\begin{array}{c}E={\displaystyle {\int }_{\Omega }}\left(K_{B}T{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{c}_i\ln c_i+\frac{1}{2}\left({\rho }_0+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{z}_{i}e{c}_i\right)\varphi \right)\text{d}x\\ +{\displaystyle \underset{i,j=1}{\overset{N}{\sum }}}{\displaystyle {\int }_{\Omega }}{\displaystyle {\int }_{\Omega }}\frac{{\epsilon }_{ij}}{2}\frac{{\left(a_i+a_j\right)}^{12}}{{\left|x-y\right|}^{12}}{c}_i\left(x\right)c_j\left(y\right)\text{d}y\text{d}x\end{array}$ (2.2)

given that $a_i$ and $a_j$ are the radii of the i^th and j^th ions and ${\epsilon }_{ij}$ becomes their coupling energy constant. The main computational challenge is the inaccuracy and inefficiency in simulation due to the effects of high frequencies in the kernel of potential term in (2.2). In order to explore the flow with high accuracy, this paper employs a band-limit function which depends on a cut-off length ( $\sigma$ ) to eliminate the high spatial frequencies and preserve the bounded spatial frequencies, see [8] [19]. The cutoff length is taken to be small parameter tending to zero for better approximation. Using this band-limit function and fourier analysis an approximate energy functional is derived which reduces the numerical complexity of the repulsive term in the LJ potential describing ion interaction.

The Energy functional in (2.2) is replaced by an equivalent approximate energy functional derived in [19] and represented as:

$\begin{array}{c}{E}_{\delta }={\displaystyle {\int }_{\Omega }}\left(K_{B}T{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{c}_i\ln c_i+\frac{1}{2}\left({\rho }_0+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{z}_{i}e{c}_i\right)\varphi \right)\text{d}x\\ +{\displaystyle \underset{i,j=1}{\overset{N}{\sum }}}\frac{{\epsilon }_{ij}}{2}{\left(a_i+a_j\right)}^{12}{S}_{\sigma }{\displaystyle \int }{c}_i\left(x\right)c_j\left(y\right)\text{d}y\text{d}x\end{array}$ (2.3)

where $S_{\sigma }=\frac{w_d}{12-d}{\sigma }^{d-12}$, the dimensional space $d\le 3$ and $w_d$ is the surface

area of the d-dimensional unit ball. Applying the energy variational method, we obtain a new mathematical model for NP equation given by:

$\frac{\partial c_i}{\partial t}=\frac{D_i}{K_{B}T}\nabla \cdot \left(c_i\nabla \frac{\delta E_{\sigma }}{\delta c_i}\right)=-\nabla \cdot J_i,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}i=1,\cdots ,N$ (2.4)

where the flux $J_i$ is given by:

$J_i=-D_i\nabla c_i-\frac{D_i{c}_i}{K_{B}T}{z}_{i}e\nabla \varphi -\frac{D_i{c}_i}{K_{B}T}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}\nabla S_{\sigma }{\epsilon }_{ij}{\left(a_i+a_j\right)}^{12}{c}_j\left(y\right)\text{d}y.$

Equation (2.4) forms a system of NP equations with steric effects, which satisfies the dissipation law given below:

$\frac{\text{d}{E}_{\sigma }}{\text{d}t}=-{\displaystyle {\int }_{\Omega }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}\frac{D_i{c}_i}{K_{B}T}{\left|\nabla \left(K_{B}T\log c_i+z_{i}e\varphi +{\mu }_i\right)\right|}^2$ (2.5)

where ${\mu }_i=\frac{\delta E_{\sigma }}{\delta c_i}$ is the chemical potential. Simplifying (2.4), we obtain rates of

concentration for the two ions $c_n$ (negative) and $c_p$ (positive) in steady state given by:

$D_n\left[\nabla \cdot \left(\nabla c_n+\frac{z_{n}e}{K_{B}T}{c}_n\nabla \varphi \right)+S_{\sigma }\nabla \cdot \left(g_{nn}{c}_n\nabla c_n+g_{np}{c}_n\nabla c_p\right)\right]=\mathrm{0,}$ (2.6a)

$D_p\left[\nabla \cdot \left(\nabla c_p+\frac{z_{p}e}{K_{B}T}{c}_p\nabla \varphi \right)+S_{\sigma }\nabla \cdot \left(g_{pp}{c}_p\nabla c_p+g_{np}{c}_p\nabla c_n\right)\right]=0.$ (2.6b)

where $g_{nn}={\epsilon }_{11}{\left(2a_{11}\right)}^{12}$, $g_{np}={\epsilon }_{12}{\left(a_1+a_2\right)}^{12}$, $g_{pp}={\epsilon }_{22}{\left(2a_{22}\right)}^{12}$. Coupling (2.6) and (2.1b) and taking ${\rho }_0=0$, we obtain the system of mPNP equations.

3. Numerical Schemes

Classical and Taylor-Hood Domain Discretisation

Case I: The study examined use Classical rectangular element models the flow where bilinear elements are prescribed on each node sharing all the three components for concentration and potential in the coupled system. Shape functions for which uses local coordinates( $\eta \mathrm{,}\xi$ ) to express these components in form of 4 noded bilinear elements $M_j$ as illustrated in Figure 1(a), representing nodes for $j=1,2,3,4$. These bilinear interpolation functions are used for both concentration, $c_n$ and $c_p$ and potential, $\varphi$ components, resulting into 12 unknown variables for each element.

$M_1=\frac{1}{4}\left[1-\xi -\eta +\xi \eta \right],-1\le \left(\xi ,\eta \right)\le 1$ (3.1)

Case II: Taylor-hood elements discretization over the domain proposes nine nodes for concentration and four nodes for potential variables where all the flow variables are defined as illustrated in Figure 1(b). Shape functions using local co-ordinates ( $\eta \mathrm{,}\xi$ ) are used to expressed these components in form of four-noded

Figure 1. The (a) Classical domain consisting of 4-nodes ( $\square$ ) providing bilinear elements for both concentration and potential components, and the (b) Taylor-hood multi-grid composed of 9-noded ( $•$ ) concentration measurements and ( $\square$ ) bilinear potential nodes.

bilinear elements given in (3.1) and nine-noded elements denoted by $N_j$ for $j=1,\cdots ,9$ an example of first node in the element is given by:

$N_1=\frac{1}{4}\xi \left(1-\xi \right)\eta \left(1-\eta \right).$ (3.2)

Thus the interpolation functions $N_j$ are used for concentration components, $c_n$ and $c_p$ while bilinear interpolation functions $M_i$ are used for potential component, $\varphi$ resulting into 20 unknown variables for each element in Taylor-hood elements.

Taylor-hood elements are established to be well-posed and numerically stable while bilinear elements are unstable [20] [21]. However, Taylor-hood imposes computational challenges due to many degrees of freedom for the variables whereas the bilinear elements easy to implement particularly when applying mesh refinement of the domain. Lastly, both elements are comparable since they produce the same order of potential interpolation.

4. Finite Element Discretization and Linearization

We employ Sobolev function to represent the domain given $C^0\left(\bar{\Omega }\right)$ to be a space of continous function at the closure of the domain, $\Omega$, $L^2\left(\Omega \right)$ a space of square integrable functions and $H_0^1\left(\Omega \right)$ being functions of $L^2\left(\Omega \right)$ whose derivatives are square integrable functions with values vanish at the boundary of the domain. These space functions are used to describe the finite dimensional

space on $\Omega \in {ℝ}^2$ to be $C_h=\left\{c_i\in C^0\left(\bar{\Omega }\right)\cap {\left(H_0^1\left(\Omega \right)\right)}^2\right\}$ and

${\Phi }_h=\left\{\Phi \in C^0\left(\bar{\Omega }\right)\cap \left(L_0^2\left(\Omega \right)\right)\right\}$, for $i=n,p$. It is therefore permissible to express parameter values $c_n$, $c_p$ and $\varphi$ are expressed in form of the basic functions of $C_h$ and ${\Phi }_h$ given by:

$c_n={\displaystyle \underset{j=1}{\overset{ndo{f}_{c_n}}{\sum }}}{N}_j{c}_{nj},c_p={\displaystyle \underset{j=1}{\overset{ndo{f}_{c_p}}{\sum }}}{N}_j{c}_{pj},\Phi ={\displaystyle \underset{j=1}{\overset{ndo{f}_{\Phi }}{\sum }}}{M}_j{\Phi }_j.$ (4.1)

where $c_{nj}$, $c_{pj}$ and ${\Phi }_j$ are the parameter measurements at the nodes of each elements and $N_j=\left\{N_1\mathrm{,}\cdots ,N_9\in C_h\right\}$ and $M_j=\left\{M_1\mathrm{,}\cdots ,M_4\in {\Phi }_h\right\}$ are the shape functions for concentration and potential respectively and whilst $ndof$ is the degree of freedom for each variable. To overcome the nonlinearity in the mPNP equations, the unknowns are written in iterative form as:

$c_n^{k+1}=c_n^k+{\hat{c}}_n,c_p^{k+1}=c_p^k+{\hat{c}}_p,{\Phi }^{k+1}={\Phi }^k+\hat{\Phi }$ (4.2)

where $c_n^k\mathrm{,}{c}_p^k$ and ${\Phi }^k$ are previous known values and ${\hat{c}}_n\mathrm{,}{\hat{c}}_p$ and $\hat{\Phi }$ are the correction values. Hence the nonlinear terms are linearized as shown in Equation (4.3):

$c_n^{k+1}\nabla c_n^k\simeq c_n^k\nabla c_n^k+{\hat{c}}_n\nabla c_n^k+c_n^k\nabla {\hat{c}}_n$ (4.3)

For novelty in computation, we let ${\alpha }_1=Z_{n}e/K_{B}T$, ${\alpha }_2=Z_{p}e/K_{B}T$, $S_1=S_{\sigma }gnn$, $S_2=S_{\sigma }gnp$ and $S_3=S_{\sigma }gpp$. Using the above notations we integrate (2.6) and (2.1b) over the domain, then by applying Gauss-Divergence theorem to the diffusive terms results into weak form of concentration and potential given by:

$\begin{array}{l}{D}_n{\displaystyle {\int }_{\Omega }}[\nabla {\hat{c}}_n+{\alpha }_1\left(c_n^k\nabla \hat{\Phi }+{\hat{c}}_n\nabla {\Phi }^k\right)+S_1\left(c_n^k\nabla {\hat{c}}_n+{\hat{C}}_n\nabla c_n^k\right)\\ +S_2\left(c_n^k\nabla {\hat{c}}_p+{\hat{c}}_n\nabla c_p^k\right)]\cdot \nabla w\text{d}\Omega \\ ={\displaystyle {\int }_{\Omega }}{f}_1\cdot w-D_n\left[\nabla c_n^k+{\alpha }_1{c}_n^k\nabla {\Phi }^k+S_1{c}_n^k\nabla c_n^k+S_2{c}_n^k\nabla c_p^k\right]\cdot \nabla w\text{d}\Omega \end{array}$ (4.4)

${\displaystyle {\int }_{\Omega }}\epsilon \nabla \hat{\Phi }\cdot \nabla v-{\displaystyle {\int }_{\Gamma }}v\cdot \frac{\partial \Phi }{\partial \zeta }-e\left({\hat{c}}_n+{\hat{c}}_p\right)\cdot v\text{d}\Omega ={\displaystyle {\int }_{\Omega }}\left(f_2-e\left(c_n^k+c_p^k\right)\right)\cdot v\text{d}\Omega$ (4.5)

given $\zeta$ is the unit outward normal and $w\in C_h\mathrm{,}v\in {\Phi }_h$ are the weight functions for concentration and potential respectively. Since achieving convergence for highly nonlinear system is computationally expensive, we refine the meshes in the domain and decouple our system to obtain convergence of the iterative numerical scheme illustrated in Figure 2 below.

5. Numerical Results and Discussion

5.1. Solutions for the PNP Problem

Case I: To validate the numerical scheme, we construct a 2D steady PNP problem with known analytical solution this is by modifying the 3D PNP equation in [22]. We start by considering analytical solution of the form (5.1a)-(5.1d) in domain, $\Omega \in {ℝ}^2$ defined as:

$\Phi \left(r\right)=a_1\cos x\cos y\mathrm{,}r\in \Omega \mathrm{,}$ (5.1a)

$c_n\left(r\right)=a_2\cos x\cos y+{\beta }_1\mathrm{,}$ (5.2b)

$c_p\left(r\right)=a_3\cos x\cos y+{\beta }_2\mathrm{,}$ (5.3c)

$F=a_4\cos x\cos y$ (5.4d)

satisfying the steady state PNP equations for two monovalent ion species $c_n$ and $c_p$ given by (5.2a)-(5.2c):

Figure 2. Flow chart for numerical simulation of mPNP equations, given $i=\left(n\mathrm{,}p\right)$.

$F_1=\nabla \cdot D_n\left(\nabla c_n+{\alpha }_1{c}_n\nabla \Phi \right)\mathrm{,}$ (5.2a)

$F_2=\nabla \cdot D_p\left(\nabla c_p+{\alpha }_2{c}_p\nabla \Phi \right)\mathrm{,}$ (5.2b)

$-\nabla \cdot \left(\epsilon \nabla \Phi \right)={\displaystyle \underset{i=1}{\overset{n\mathrm{,}p}{\sum }}}{z}_{i}e{c}_i\mathrm{.}$ (5.2c)

when taking $D_n=1$, $D_p=0.8$, ${\alpha }_1=-1$, ${\alpha }_2=1$, $\left(a_1,a_2,a_3,a_4\right)=\left(0.2,0.2,0.4,80\right)$ and $\left({\beta }_1,{\beta }_2\right)=\left(0.3,-0.3\right)$. Given $\epsilon ={\epsilon }_r{\epsilon }_0$, ${\epsilon }_0$ is the vacuum permittivity and ${\epsilon }_r$ is the relative permittivity. Prescribing Neumann boundary condition for concentration in two opposite vertical walls of the channel and Dirichlet boundary conditions for exit and entry of the channel $\Omega =\left[0,1\times 0,1\right]$ walls satisfying:

$\left(\nabla c_n-{\alpha }_1{c}_n\nabla \Phi \right)\cdot \nu =0,$ (5.3a)

$\left(\nabla c_p+{\alpha }_1{c}_n\nabla \Phi \right)\cdot \nu =0,$ (5.3b)

$\epsilon \nabla \Phi \cdot \nu =0,$ (5.3c)

$c_i\left(r\right)=k_i\left(r\right),{\Phi }_r=h\left(r\right),r\in \partial \Omega .$ (5.3d)

where $\nu$ is the unit outward normal.

The numerical solution of PNP equations (5.2a)-(5.2c) was in agreement with the analytical solution for Concentration and Potential as illustrated in Figure 3 and Figure 4. Increasing the number of elements, N improved the accuracy of the numerical solution such that when $N=16\times 16$, the approximate solution were established to be very accurate hence validating the numerical scheme.

5.2. Solutions for mPNP Problem

In this subsection, we consider a two dimensional steady state mPNP equations given by (2.6) and (2.1b) for two ion species denoted by n and p, with valencies taken to be $z_n=-1$ and $z_p=1$, and radii $a_1=1.5$ Å and $a_2=2.0$ Å, respectively. The constant diffusion coefficients was taken as $D_n=2.0305\times {10}^{-5}$ and $D_p=1.335\times {10}^{-5}$ for negative and positive ions, respectively. The flow is assumed to be in the y direction at $x=0.5$. The impact of steric components $S_1\mathrm{,}{S}_2$ and $S_3$ on concentration of two monovalent ion species is determined and discussed for Taylor-hood and bilinear elements.

Dirichlet conditions prescribed at the channel walls for concentration and potential in the domain, $\Omega =\left[\mathrm{0,1}\times \mathrm{0,1}\right]$, no flux boundary conditions are imposed for both concentration and potential, respectively and electro-neutrality condition is assumed for the charge densities given by:

$\begin{array}{l}{c}_n\left(x,0\right)=c_p\left(x,0\right)=g1\hfill \\ c_n\left(x,1\right)=c_p\left(x,1\right)=0\hfill \\ \Phi \left(x,0\right)=g2,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\Phi \left(x,1\right)=0\hfill \\ J_i\cdot \nu =0,\left(\epsilon \nabla \Phi \right)\cdot \nu =0\hfill \\ z_n{c}_n+z_p{c}_p=0\hfill \end{array}\}\text{on}\partial \Omega ,i=n,p$ (5.4)

where $g1>0,g2>0$ are constants, $J_i$ is the concentration flux and $\nu$ is the

Figure 3. Contours for (a) analytical and (b) numerical Concentration, c_n satisfying the PNP equations (5.2a)-(5.2c).

unit outward normal.

Case II: A four-noded rectangular element, each sharing nodal points for potential and concentration components was utilized to obtain the following results. The three steric effects $S_1\mathrm{,}{S}_2$ and $S_3$ are varied and analyzed for concentration profile with the aim to demonstrate their effect on the positive and negative ion flow in an electro-neutral biological channel environment under the

Figure 4. Potential, $\Phi$ contours for (a) analytical and (b) numerical solutions of the PNP equations.

boundary condition in (5.4).

Results illustrated in Figure 5 demonstrate effect of repulsive forces of cations in the flow whence other steric components are held constant. On the other hand, Figure 6 and Figure 7 accounts for the attractive force of the constituent

Figure 5. Variation of ionic concentration c_n and c_p with finite size effects, S₃ and PNP equation across the domain, at $x=0.5$.

opposite ions and repulsive of anions, respectively, all of which are responsible for ion flow fluctuations.

Figure 6. Variation of ionic concentration with finite size effects, S₂ across the domain at $x=0.5$.

The above three critical forms in variation implies that radius size of an ion plays a role in ion interaction and competition as in the subsequent analysis.

Figure 7. Variation of ionic concentration with finite size effects, S₁ across the domain at $x=0.5$.

Upon increasing the positive steric effects S₃ the repulsive forces between ion increases resulting into reduction in cations being selected to through the channel see, Figure 6(b). These repulsive forces are very strong hence prevents of anions as displayed in Figure 6(b). Varying the attractive element of the steric effect, S₂ increases the attractive forces between the ions resulting in gradual increase in flow and accumulation of anions and accelerated decrease of the cations in the channel as shown in Figure 7. This may be due to frequent collision between cations which have larger radius hence reaches equilibrium faster than anions. Lastly, insignificant variation of the negative steric effects S₁ increases the repulsive forces between the anions enhancing steady reduction in their diffusion and concentration in the channel, see Figure 8(a) whereas the flow of cations remain constant, as evident in Figure 8(b). Consequently, we realize significant contribution of both attractive and repulsive element of steric effect compared to ion concentration. But overally, contribution of the increase in the attractive component of the steric effect S₂ in the flow pronounces anion and diminishes cation flows significantly as in Figure 7(a) and Figure 7(b) respectively.

It can, therefore, be deduced that the repulsive components plays the least role in the selectivity of ion species while the attractive components have more impact on the flow of ions. However flow of the anions is enhanced in all the situations compared to the cations upon varying the repulsive components S₃ and S₁. This may be as a result of the size/radius of the ion. The diffusion of cations is also observed to be faster than anions in all the variations, this may be as a result of bulk mobility of the anions which is higher than the cations in the channel.

Case III: Taylor-hood multi-grid with 9 nodes rectangular elements for each concentration components and 4 noded rectangular elements for potential components was adopted to discretize and solve the mPNP system of equations. Given, specific values of steric effects S₁, S₂ and S₃, the estimated concentration profiles c_n and c_p obtained are similar and showed convergence for gradual increase in the number of iterations as illustrated in Table 1. Whereas execution of infinitesimal variation of steric effect, S₁ results into insignificant variation in the approximate solutions of the flow components, see Figure 8. This loss of accuracy and convergence is largely because of the instability in the system due to nonlinearity and loss of information from non shared nodes.

The results obtained with mPNP equations demonstrated the impact of variation of finite size on ion interaction in the channel, which is not viable with PNP equations. Simulation with bilinear elements exhibited numerically stable solutions of the mPNP model compared to Taylor-hood elements, this may be as a consequence of few degrees of freedom of potential nodes updating the concentration nodes in each element in the domain. Further, results with Taylor-hood elements introduced oscillations in the approximate solution and some negative concentration with varied steric component. Therefore, it is importance to address this challenge in order to obtain desired results. Consequently, we investigated the convergence of the iterative scheme after number of iterations, as illustrated in Table 1 and results with Taylor-hood elements displayed instability in convergence compared to bilinear elements and PNP equations.

Figure 8. Concentration profile for (a) c_n and (b) c_p for various, S₁ and PNP equations when $x=0.5$.

Table 1. Convergence of exact PNP, bilinear (cFEM) and Taylor-hood (TH) elements after iterations I.

6. Conclusion

The main objective of the study was to develop PNP with steric effects consisting of LJ hard-sphere potential which was modified using a band limit function to reduce the complexity of computation. A two-dimensional steady-state numerical solution of the mPNP system of equations showing the effects of variation of steric effects on ion flow and concentration is discussed for a rectangular channel. In fact, we have observed the effects of repulsive and attractive steric forces on ion flow and deduced the role played by radial size in selectivity. The Mixed Finite element approach based on Taylor-hood and bilinear elements enabled the establishment of distinction between the flows in relation to steric components qualitatively. We validated our numerical scheme by obtaining analytical and numerical solutions on the 2D PNP equation and convergence of the scheme compared to mPNP equations. Lastly, it is fundamental to note that an anion has a smaller size compared to a cation, therefore, easing selectivity and flow in the channel as established in the study. Effects of potential variation are also worth examining for such a flow in addition to computational efficiency in triangular elements in Taylor-hood method.

Conflicts of Interest

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

References

[1]	Chen, D. and Wei, G.W. (2017) A Review of Mathematical Modeling, Simulation and Analysis of Membrane Channel Charge Transport. In: Reference Module in Life Sciences, Elsevier, Amsterdam, 1-42. https://doi.org/10.1016/B978-0-12-809633-8.12044-8
[2]	Zheng, Q. and Wei, G.W. (2011) Poisson-Boltzmann-Nernst Planck Model. The Journal of Chemical Physics, 134, Article ID: 194101. https://doi.org/10.1063/1.3581031
[3]	Eisenberg, B. (1998) Ionic Channels in Biological Membranes-Electrostatic Analysis of Natural Nanotube. Contemporary Physics, 39, 447-466. https://doi.org/10.1080/001075198181775
[4]	Yang, Z.C., Van, T.A., Straaten, D., Ravaioli, U. and Liu, Y. (2005) A Coupled 3-D PNP/ECP Model for Ion Transport in Biological Ion Channels. Journal of Computational Electronics, 4, 167-170. https://doi.org/10.1007/s10825-005-7131-8
[5]	Flavell, A., Machen, M., Eisenberg, B., Kabre, J., Liu, C. and Li, X.F. (2013) A Conservative Finite Difference Scheme for Poisson Nernst-Planck Equations. Journal of Computational Electronics, 13, 235-249. https://doi.org/10.1007/s10825-013-0506-3
[6]	Meng, D., Zheng, B., Lin, G. and Sushko, M.L. (2014) Numerical Solution of 3D Poisson Nernst-Planck Equations Coupled with Classical Density Function Theory for Modeling Ion and Electron Transport in a Confined Environment. Communications in Computational Physics, 16, 1298-1322. https://doi.org/10.4208/cicp.040913.120514a
[7]	Siddiqua, F., Wang, Z.M. and Zhou, S.G. (2017) A Modified Poisson-Nernst-Planck Model with Excluded Volume Effects: Theory and Numerical Implementation. Communications in Mathematical Sciences, 16, 26 p.
[8]	Horng, T.L., Lin, T.C., Liu, C. and Eisenberg, B. (2012) PNP Equations with Steric Effects: A Model of Ion Flow through Channels. The Journal of Physical Chemistry B, 116, 11422-11441.
[9]	Lu, B.Z. and Zhou, Y.C. (2011) Poisson Nernst Equations for Simulating Biomolecular Diffusion-Reaction Processes II: Size Effects on Ionic Distributions and Diffusion-Reaction Rates. Biophysical Journal, 100, 2475-2485. https://doi.org/10.1016/j.bpj.2011.03.059
[10]	Hyon, Y., 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
[11]	Qian, Y., Liu, X.J., Chen, M.X. and Lu, B.Z. (2016) A Local Approximation of Fundamental Measure Theory Incorporated into Three Dimensional Poisson-Nernst Planck Equations to Account for Hard Sphere Repulsions among Ions. Journal of Statistical Physics, 163, 156-174. https://doi.org/10.1007/s10955-016-1470-7
[12]	Vladmir, Y., Kiselev, M., Lobanov, I., Marenduzzo, D. and Goryachev, B. (2011) Lateral Dynamics of Charged Lipids and Peripheral Proteins in Spatially Heterogeneous Membranes: Comparison of Continuous and Monte Carlo Approaches. Journal of Chemical Physics, 135, Article ID: 155103. https://doi.org/10.1063/1.3652958
[13]	Gao, H.D. and Sun, P.T. (2018) A Linearized Local Conservative Mixed Finite Element Method for Poisson-Nernst-Planck Equations. Journal of Scientific Computing, 77, 793-817.
[14]	Burger, M., Schlake, B. and Wolfram, M.T. (2012) Nonlinear Poisson-Nernst-Planck Equations for Ion Flux through Confined Geometries. Nonlinearity, 25, 961-990. https://doi.org/10.1088/0951-7715/25/4/961
[15]	Gwecho, A.M., Wang, S. and Mboya, O.T. (2020) Existence of Approximate Solution for Modified Poisson Nernst-Planck Equation Describing Ion Flow in Cell Membranes. American Journal of Computational Mathematics, 10, 473-484. https://doi.org/10.4236/ajcm.2020.103027
[16]	Abidha, M., Wang, S. and Onyango, T. (2020) Numerical Approach for Determining Impact of Steric Effects in Biological Ion Channel. Journal of Computational and Applied Mathematics, 9, 1-4.
[17]	Eisenberg, B. (2011) Mass Action in Ionic Solutions. Chemical Physics Letters, 511, 1-6. https://doi.org/10.1016/j.cplett.2011.05.037
[18]	Tu, B., Xie, Y., Zhang, L. and Lu, B. (2015) Stabilized Finite Element Methods to Simulate the Conductance of Ion Channels. Computer Physics Communications, 188, 131-139. https://doi.org/10.1016/j.cpc.2014.11.018
[19]	Lin, T.C. and Einsenberg, B. (2013) A New Approach to the Lennard-Jones Potential and a New Model: PNP-Steric Equations. Communications in Mathematical Sciences, 12, 149-173.
[20]	Brezzi, F. and Fortin, M. (1991) Mixed and Hybrid Finite Element Methods. Springer, New York. https://doi.org/10.1007/978-1-4612-3172-1
[21]	White, J.A. and Borja, R.I. (2008) Stabilized Low-Order Finite Elements for Coupled Solid-Deformation/Fluid-Diffusion and Their Application to Fault Zone Transients. Computer Methods in Applied Mechanics and Engineering, 197, 4353-4366. https://doi.org/10.1016/j.cma.2008.05.015
[22]	Zheng, Q., Chen, D. and Wei, G.W. (2011) Second Order Poisson-Nernst-Planck Solver for Ion Channel Transport. Journal of Computational Physics, 230, 5239-5262. https://doi.org/10.1016/j.jcp.2011.03.020