A free energy satisfying finite difference method for Poisson–Nernst–Planck equations
Abstract.
In this work we design and analyze a free energy satisfying finite difference method for solving PoissonNernstPlanck equations in a bounded domain. The algorithm is of second order in space, with numerical solutions satisfying all three desired properties: i) mass conservation, ii) positivity preserving, and iii) free energy satisfying in the sense that these schemes satisfy a discrete free energy dissipation inequality. These ensure that the computed solution is a probability density, and the schemes are energy stable and preserve the equilibrium solutions. Both one and twodimensional numerical results are provided to demonstrate the good qualities of the algorithm, as well as effects of relative size of the data given.
Key words and phrases:
Poisson equation, NernstPlanck equation; free energy; positivity1991 Mathematics Subject Classification:
35K20, 65M06, 65M12, 82C31.1. Introduction
In this paper, we are interested in constructing a free energy satisfying numerical method for solving the initial boundary value problem for the Poisson–Nernst–Planck (PNP) equations,
(1.1) 
where is the concentration of ion species, denotes a connected closed domain with smooth boundary , is the potential governed by the Poisson equation which is necessary to determine the electrostatic field, and n is the unit outward normal vector. Subject to the given initial and boundary conditions, the compatibility condition
(1.2) 
is necessarily to be imposed for solvability of the problem. By freeenergy satisfying we mean that the free energy dissipation law is satisfied at the discrete level.
The PNP equations describe the diffusion of ions under the effect of an electric field that is itself caused by those same ions. The system couples the NernstPlanck (NP) equation (which describes the drift of ions in a potential gradient by Ohm’s law and diffusion of ions in a concentration gradient by Fick’s law) and the Poisson equation (which relates charge density with electric potential). This system of equations for multiple species has been extensively used in the modeling of semiconductors (see e.g., [32]), and the membrane transport in biological ion channels (see e.g., [14]).
The PNP system can hardly be solved analytically. The main difficulty arises from the nonlinear coupling of the electrostatic potential and concentrations of chemical species. When the physical domain has a simple geometry, a semiexplicit formula was derived in [26] for the steadystate solution; the existence and stability of the steadystate solution was established a long while ago [21] in the study of the steady Van Roostbroeck model in semiconductors. It has been proved by H. Gajewski and K. Gärtner in [18] that the solution to the driftdiffusion system converges to the thermal equilibrium state as time becomes large if the boundary conditions are in thermal equilibrium. The keypoint of the proof is an energy estimate with the control of the free energy dissipation. Long time behavior was also studied in [5], and further in [1, 4] with refined convergence rates.
In the past decade a growing interest in PNP systems has been driven mainly by experimental and numerical advances. Computational algorithms have been constructed for both simple onedimensional settings and complex threedimensional models in various chemical/biological applications, and have been combined with the Brownian Dynamics simulations; cf. [13, 22, 9, 20, 17, 11, 30, 26, 3, 33, 34, 29, 36, 37, 39]. Many of these existing algorithms are introduced to handle specific settings in complex applications, in which one may encounter different numerical obstacles, such as discontinuous coefficients, singular charges, geometric singularities, and nonlinear couplings to accommodate various phenomena exhibited by biological ion channels [40]. For instance, [36, 37] developed a finite difference scheme for liquid junction and ion selective membrane potentials, in which authors used a fully implicit discretization scheme and the NewtonRaphson solver for the resulting linear system. In [30], a 3D finite element method for the system with a singular (pointlike) charge outside of the concentration region was developed. A second order PNP solver was developed in [39] in a realistic ionchannel context with the Dirichlet boundary condition. In spite of many existing computational studies, rigorous numerical analysis seems to be still lacking.
Our objective is to construct and analyze an explicit secondorder PNP algorithm to incorporate main mathematical features of the PNP system so that the numerical solution remains faithful for long time simulations, i.e., the numerical solution possesses desired properties including conservation of ions, positivity of concentration and dissipation of the free energy. Therefore we consider only the standard PNP equations of form (1.1), while the reader may find intensive discussions of the modeling aspect of the PNP system in the literature.
The main properties of the solution to (1.1) are the nonnegativity principle, the mass conservation and the free energy dissipation, i.e.,
(1.3)  
(1.4)  
(1.5) 
where the free energy is defined by
If does not depend on time, we can use a modified functional so that
These properties are also naturally desired for numerical methods solving (1.1). In this paper, we develop such a method. We will demonstrate that these properties of the numerical methods could be critical in obtaining the longtime behavior of the solutions.
It is difficult for numerical schemes to preserve all three properties for PNP equations exactly at the discrete level. A recent effort toward this direction is found in [16], where the authors present an implicit second order finite difference scheme with a simple iteration so that the energy dissipation law is approximated closely.
The free energy dissipation in time is also the driving force so that the large time behavior of the solution to (1.1) is governed by the steadystate solution. In fact, the steadystate solution is necessarily of the form , which gives , and can be further determined by integration of the Poisson equation while using the Neumann boundary data. In other words, the PNP system near steady states is close to a Poisson–Boltzmann equation (PBE)
Therefore, the numerical method presented in this paper may be used as an iterative algorithm to numerically compute such a nonlocal PBE. We refer to [23] for a rigorous study of a nonlocal Poisson–Boltzmann type equation with a small dielectric parameter as derived from the twospecies PNP system. In general, the PBE describes the electrostatic interaction and ionic density distributions of a solvated system at the equilibrium state. Due to the effectiveness of the PBE for applications in chemistry and biophysics, a large amount of literatures and many solution techniques have been produced in this area and directed to studies of diverse biological processes, see e.g. [31, 10, 12, 15, 19, 24, 25, 35].
1.1. Related models
In a larger context, the concentration equation falls into the general class of aggregation equations with diffusion
(1.6) 
where . Such a model has been widely studied in applications such as biological swarms [7, 8, 38] and chemotaxis [2, 6]. For chemotaxis, a wide literature exists in relation to the KellerSegel model (see [2, 6] and references therein). The lefthandside in (1.6) represents the active transport of the density associated to a nonlocal velocity field . The potential is assumed to incorporate attractive interactions among individuals of the group, while repulsive (anticrowding) interactions are accounted for by the diffusion in the righthandside.
Of central role in studies of model (1.6), and also particularly relevant to the present research, is the gradient flow formulation of the equation with respect to the free energy
(1.7) 
The free energy contains both entropic part and the interaction part. There is a vast literature on entropic schemes for kinetic equations such as FokkerPlanck type equations. For these equations, information carried by the probability density becomes less and less as time evolves, the probability density is expected to converge to the equilibrium solution in a closed system regardless of how initial data are distributed. The entropy dissipation in time is the underlying mechanism for this phenomenon. The entropy satisfying methods recently developed in [27, 28] are the main motivation for this paper.
1.2. Contents
This paper is organized as follows: in Section 2, we describe the algorithm for the one dimensional case. Theoretical analysis for both semidiscrete and full discrete schemes is provided. We also discuss how the Poisson equation is solved. Numerical results of both one and two dimensions are presented in Section 3. Finally, in Section 4, concluding remarks are given.
2. The numerical method
2.1. Reformulation
Following [27], we formally reformulate the system by setting
to obtain
(2.8)  
(2.9) 
With this reformulation, an iterative algorithm may be designed as follows: given we obtain by solving the Poisson equation
which determines ; we further update by solving (2.8), (2.9).
We shall describe our numerical algorithm in one dimensional setting only, with one or multispecies. The algorithm is extensible to high dimensional case in a straightforward dimension by dimension manner.
In the one dimensional case with , the above two steps correspond to solving the following two sets of problems:
(2.10)  
(2.11)  
(2.12)  
(2.13)  
(2.14) 
Note that we set and to different constants. Integrating the Poisson equation and using the mass conservation for we obtain
which is exactly the compatibility condition stated in (1.2).
If are independent of time, we set
A direct calculation gives
We now describe our algorithm by first partitioning the domain with , , and interior grid points
2.2. Algorithm

We use to approximate and to approximate , respectively. Given , we compute the potential by
(2.15) where , and . For definiteness, we set at any time to single out a particular solution since is unique up to an additive constant.

With the above obtained , the semidiscrete approximation of the concentration satisfies
(2.16) and , and .

Discretize uniformly: , so and satisfy
(2.17)
For the system case with the vector unknown , the PNP system in a simple form is
(2.18a)  
(2.18b) 
where is the charge of . The onedimensional algorithm needs to be modified: (2.15) is replaced by
(2.19) 
where , in order to solve the Poisson equation (2.18b), and each is obtained from solving
(2.20) 
where in the formula for the potential is replaced by . For the multidimensional case, the algorithm can be applied in a dimension by dimension manner. We illustrate these cases by numerical examples later.
Remark 2.1.
The above algorithm can be easily modified to solve the PNP equations with more physical parameters:
(2.21a)  
(2.21b) 
where is the ion density for the th species, is the charge of , is the diffusion constant, is the Boltzmann constant, is the absolute temperature, is the permittivity, is the electrostatic potential, is the permanent (fixed) charge density of the system, and is the number of ion species.
Remark 2.2.
It is also possible to generalize the schemes presented in this work to second order finite difference schemes on nonuniform meshes, yet the analysis would appear more complicated.
2.3. Properties of the numerical method
The numerical solution obtained from the above algorithm has some desired properties as stated in the following.
Theorem 2.1.

The discrete concentration remains positive in time: if , then
provided the condition where
(2.24) 
The semidiscrete free energy
satisfies
(2.25) therefore nonincreasing.
Proof.

A direct calculation using gives
We proceed with
(2.28) since and for any and . Next we use the boundary conditions and to obtain
(2.29) Putting all together leads to (2.25).
∎
Remark 2.3.
For multispecies case governed by (2.18), the algorithm still preserves three desired properties. Mass conservation follows from (2.20) for each so that
For the positivity property, we need to bound , which in multispecies case becomes
where . Using the fact that we have
These bounds of ensure that provided that
where is given in (2.24). This sufficient condition is consistent with the onespecie case . For the free energy dissipation property, we consider the semidiscrete functional of the form
Following the calculation in (2.28) for each , and (2.29) using (2.19) for the potential part, we arrive at
(2.30) 
where .
2.4. On solving the Poisson equation
One of the main numerical tasks is to solve the Poisson equation
in the first step of our algorithm, subject to the Neumann conditions. In one dimensional case, we use the numerical boundary condition and to construct the tridiagonal matrix and the right hand side source term
so that the corresponding linear system becomes with . However, the above defined is singular since the solution to the Poisson equation with Neumann boundary conditions is unique up to an additive constant. We pick one definite solution by setting , i.e., we set
In this new formulation, the surface charge is implicitly used because of the compatible condition
which can be seen from the sum of (2.15) over . In the old formulation, the first equation gives after we set . In the new formulation, the same can be obtained by summing all equations except the first one:
In two dimensions of a rectangle domain with a uniform partition of and , we first arrange and into an array row by row, i.e.,
where and approximate and , respectively.
With the Neumann boundary condition , and in a square domain(), we solve the linear system , where


and
Note that if , the matrices , , , and are no longer squares, and need to be modified accordingly. Note that the matrix is a constant matrix in our setting.
3. Numerical examples
We now test our method in various settings.
3.1. One dimensional numerical tests
Single species
We first consider a neutral system with , , and , which satisfies the compatibility condition (1.2). The computational domain is and initial conditions are
Table 3.1 shows the numerical convergence at time . The numerical solutions with are taken to be discrete reference solutions, and a piecewise cubic spline interpolation is used to obtain the reference solution and . The error of the numerical solution with step size is then defined by , where are the grid points. The convergence orders are defined by . The same definition applies to . We observe that the spacial accuracy of the numerical scheme is of second order in both and for all three initial data, including the discontinuous data.
Fig. 3.1 shows the initials and final steadystate solutions, where we notice that solutions with all initials converge to the same steadystate solution given the same stop criteria, if the boundary data and the total density are the same. The initial in (3) is set to be zero in , the numerical solution will maintain positivity all the time and handle the discontinuity at well. Fig. 3.2 shows the free energy decay in time, where the free energy is truncated in for better comparisons. By adopting the same convergence criteria in the energy , we observe that it takes the longest time for the initial data in (3) to evolve into the steady state. This is as anticipated, since the initial in (3) is farther away from the steady state than other initial data. Eventually the free energy converges to for all three cases. Long time solution behavior is also tested, and our results show that the energy stays at even at for all three cases.
h  error in  order  error in  order 
0.2  0.0026321    0.00041461   
0.1  0.00065636  2.0037  0.00012431  1.7378 
0.05  0.0001639  2.0017  3.4138e005  1.8645 
0.025  4.2017e005  1.9637  9.1157e006  1.9049 
0.0125  1.0026e005  2.0672  2.224e006  2.0352 
0.2  0.002514    0.00030611   
0.1  0.00060573  2.0532  0.0001046  1.5492 
0.05  0.00015138  2.0005  2.9183e005  1.8416 
0.025  4.1626e005  1.8627  6.7804e006  2.1057 
0.0125  9.9647e006  2.0626  1.6558e006  2.0338 
0.00625  2.0458e006  2.2842  3.277e007  2.3371 
0.1  0.00082461    0.00014438   
0.05  0.00020753  1.9904  3.9146e005  1.8829 
0.025  6.347e005  1.7092  1.2189e005  1.6833 
0.0125  1.5296e005  2.053  2.978e006  2.0331 
0.00625  3.2147e006  2.2504  6.2711e007  2.2476 
Multiple species
We consider a neutral system with two species on , e.g.,
In a domain [0,1], we pick and and , , and in order to satisfy compatibility condition (1.2). Initial conditions are
Fig. 3.3 shows the initial conditions and their steadystate solutions. We also observe that in all cases the numerical solutions converge to the same steadystate solution for the above fixed , and , and the free free energy converges to .
3.2. Two dimensional numerical tests
Single species
We now test our algorithm in a two dimensional setting. Let be the unit square centered at . The initial and boundary conditions are
Note that the compatibility condition (1.2) is satisfied for both cases.
Table 3.2 shows the numerical error and the convergence order for the initial case (1), where a reference solution is computed at . We observe the convergence rate at order two for and slightly worse for . Fig. 3.4 shows the steadystate solution of different initial conditions. For the case 1, we also display the evolution of the concentration (from bottom to top as time evolves). The evolution of the potential is not shown since it does not change too much visually.
h 
error in  order  error in  order 

0.2000  0.034449    0.0005512   
0.1000  0.009046  1.9291  0.00019733  1.482 
0.0500  0.0022608  2.0004  6.7967e005  1.5377 
0.0250  0.00053914  2.0681  1.867e005  1.8641 
0.0125  0.00010783  2.3219  5.9918e006  1.6397 
4. Concluding remarks
In this paper, we have investigated the PoissonNernstPlanck equation which is of mean field type model for concentrations of chemical species, with our focus on the development of a free energy satisfying numerical method for the PNP equation subject to zero flux for the chemical concentration and nontrivial flux for the potential on the boundary. We constructed simple, easytoimplement conservative schemes which preserve equilibrium solutions, and proved that they satisfy all three desired properties of the chemical concentration: the total concentrations is conserved exactly, positivity of the chemical concentrations is preserved under a mild CFL condition , and the free energy dissipation law is satisfied exactly at the semidiscrete level. The analysis is for one and multiple ionic species and for onedimensional problems. But numerical tests are given for both onedimensional and twodimensional systems, as well as for multiple species.
Although this work makes good progress in constructing and analyzing an accurate method for solving the PoissonNernstPlanck equations numerically, there remain many challenges, which we wish to address in future work. First, we would like to extend the present numerical method and analytical results herein to a higher order discontinuous Galerkin method. Second, for most ion channels, the appropriate boundary conditions are those with the nonzero flux across the boundary. In such settings, mass conservation is no longer valid, we will investigate the possibility to preserve the solution positivity and the energy dissipation property at the discrete level.
Acknowledgments
The authors thank Bo Li for stimulating discussions of the PNP system. Liu’s research was partially supported by the National Science Foundation under grant DMS 1312636.
Footnotes
 The same h is used for both x and ydirections.
References
 A. Arnold, P. Markowich and G. Toscani. On large time asymptotics for driftdiffusionPoisson systems. Transport Theory and Statistical Physics, 29(35):57158, 2000.
 J. Bedrossian, N. Rodríguez and A. L. Bertozzi. Local and global wellposedness for aggregation equations and PatlakKellerSegel models with degenerate diffusion. Nonlinearity, 24:168–714, 2011.
 S. Berneche and B. Roux. A microscopic view of ion conduction through the K+ channel. Proc. Natl. Acad. Sci. U.S.A., 100:8644, 2003.
 P. Biler and J. Dolbeault. Long time behavior of solutions to NernstPlanck and DebyeHückel driftdiffusion systems. Annales Henri Poincaré, 1(3):461–472, 2000.
 P. Biler, W. Hebisch and T. Nadzieja. The Debye system: existence and large time behavior of solutions. Nonlinear Anal., 23:118–209, 1994
 A. Blanchet, J. A. Carrillo and P. Laurencot. Critical mass for a PatlakKellerSegel model with degenerate diffusion in higher dimensions. Calc. Var. Partial Differential Equations, 35(2):133–168, 2009.
 M. Burger, V. Capasso and D. Morale. On an aggregation model with long and short range interactions. Nonlinear Analysis: Real World Applications, 8:939–958, 2007.
 M. Burger and M. D. Francesco. Large time behavior of nonlocal aggregation models with nonlinear diffusion. Netw. Heterog. Media, 3(4):749–785, 2008.
 A. E. Cardenas, R. D. Coalson and M. G. Kurnikova. Threedimensional PoissonNernstPlanck theory studies: influence of membrane electrostatics on gramicidin A channel conductance. Biophys. J., 79:80, 2000.
 J. Che, J. Dzubiella, B. Li and J. A. McCammon. Electrostatic free energy and its variations in implicit solvent models. J. Phys. Chem. B, 112:3058–3069, 2008.
 B. Corry, S. Kuyucak and S. H. Chung. Tests of continuum theories as models of ion channels. ii. PoissonNernstPlanck theory versus brownian dynamics. Biophys. J., 78:2364–2381, 2000.
 M. E. Davis and J. A. McCammon. Electrostatics in biomolecular structure and dynamics. Chem. Rev., 90:509–521, 1990.
 R. Eisenberg. Computing the field in proteins and channels. J. Membr. Biol., 150:1–25, 1996.
 R. Eisenberg. Ion channels in biological membranes: electrostatic analysis of a natural nanotube. Contemp. Phys., 39:447–, 1998.
 M. Fixman. The PoissonBoltzmann equation and its application to polyelecrolytes. J. Chem. Phys., 70:4995–5005, 1979.
 A. Flavell, M. Machen, R. Eisenberg, J. Kabre, C. Liu and X. Li. A conservative finite difference scheme for PoissonNernstPlanck equations. J. Comput. Electron., 15:1–15, 2013
 S. Furini, F. Zerbetto and S. Cavalcanti. Application of the PoissonNernstPlanck theory with spacedependent diffusion coefficients to KcsA. Biophys. J., 91:3162–3169, 2006.
 H. Gajewski and K. Gärtner. On the discretization of Van Roosbroeck’s equations with magnetic field. Z. Angew. Math. Mech., 76(5):247–264, 1996.
 P. Grochowski and J. Trylska. Continuum molecular electrostatics, salt effects and counterion binding—A review of the PoissonBoltzmann model and its modifications. Biopolymers, 89:93–113, 2008.
 U. Hollerbach, D. P. Chen, D. D. Busath and R. Eisenberg. Predicting function from structure using the PoissonNernstPlanck equations: sodium current in the gramicidin A channel. Langmuir, 16:5509–5514, 2000.
 J. W. Jerome. Consistency of semiconductor modeling: an existence/stability analysis for the stationary van Boosbroeck system. SIAM J. Appl. Math., 45:565–590, 1985.
 M. G. Kurnikova, R. D. Coalson, P. Graf and A. Nitzan. A lattice relaxation algorithm for threedimensional PoissonNernstPlanck theory with application to ion transport through the gramicidin A channel. Biophys. J., 76:642–656, 1999.
 C. Lee, H. Lee, Y. Hyon, T. Lin and C. Liu. New PoissonBoltzmann type equations: onedimensional solutions. Nonlinearity, 22(2):431–458, 2011.
 B. Li. Continuum electrostatics for ionic solutions with nonuniform ionic sizes. Nonlinearity, 22:811–833, 2009.
 B. Li. Minimization of electrostatic free energy and the PoissonBoltzmann equation for molecular solvation with implicit solvent. SIAM J. Math. Anal., 40(6):2536–2566, 2009.
 B. Li, B. Lu, Z. Wang and J.A. McCammon. Solutions to a reduced PoissonNernstPlanck system and determination of reaction rates. Physica A: Statistical Mechanics and its Applications, 389(7):1329–1345, 2010.
 H. Liu and H. Yu. An entropy satisfying conservative method for the FokkerPlanck equation of the finitely extensible nonlinear elastic dumbbell model. SIAM J. Numer. Anal., 50(3):1207–1239, 2012.
 H. Liu and H. Yu. The entropy satisfying discontinuous Galerkin method for FokkerPlanck equations, with applications to the finitely extensible nonlinear elastic dumbbell model. SIAM J. Numer. Anal., in review, 2013.
 C. L. Lopreore, T. M. Bartol, J. S. Coggan, D. X. Keller, G. E. Sosinsky, M. H. Ellisman and T. J. Sejnowski. Computational modeling of threedimensional electrodiffusion in biological systems: Application to the node of Ranvier. Biophys. J., 95:2624–2635, 2008.
 B. Z. Lu, M. J. Holst, J. A. McCammon and Y. C. Zhou. PoissonNernstPlanck equations for simulating biomolecular diffusionreaction processes I: Finite element solutions. J. Comput. Phys., 229(19):6979–6994, 2010.
 B. Z. Lu, Y. C. Zhou, G. A. Huber, S. D. Bond, M. J. Holst and J. A. McCammon. Electrodiffusion: A continuum modeling framework for biomolecular systems with realistic spatiotemporal resolution. J. Chem. Phys., 127:135102, 2007.
 P.A Markowich, C.A. Ringhofer and C. Schmeiser. Semiconductor Equations. Springer, New York, 1990.
 S. Y. Noskov, W. Im and B. Roux. Ion permeation through the hemolysin channel: Theoretical studies based on brownian dynamics and PoissonNernstPlank electrodiffusion theory. Biophys. J., 87:2299, 2004.
 B. Roux, T. Allen, S. Berneche and W. Im. Theoretical and computational models of biological ion channels. Q. Rev. Biophys., 37:15, 2004.
 K. A. Sharp and B. Honig. Calculating total electrostatic energies with the nonlinear PoissonBoltzmann equation. J. Phys. Chem., 94:7684–7692, 1990.
 T. Sokalski and A. Lewenstam. Application of NernstPlanck and Poisson equations for interpretation of liquidjunction and membrane potentials in realtime and space domains. Electrochemistry Communications, 3(3):107–112, 2001.
 T. Sokalski, P. Lingenfelter and A. Lewenstam. Numerical solution of the coupled NernstPlanck and Poisson equations for liquid junction and ion selective membrane potentials. The Journal of Physical Chemistry B, 107(11):2443–2452, 2003.
 C. M. Topaz, A. L. Bertozzi and M. A. Lewis. A nonlocal continuum model for biological aggregation. Bull. Math. Bio., 68:1601–623, 2006.
 Q. Zheng, D. Chen and G.W. Wei. Secondorder PoissonNernstPlanck solver for ion channel transport. J Comput Phys., 230(13):5239–5262, 2011.
 G.W. Wei, Q. Zheng, Z. Chen and K. Xia. Variational multiscale models for charge transport. SIAM Rev., 54:699–754, 2012.