A free energy satisfying finite difference method for Poisson–Nernst–Planck equations
In this work we design and analyze a free energy satisfying finite difference method for solving Poisson-Nernst-Planck 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 two-dimensional 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, Nernst-Planck equation; free energy; positivity
1991 Mathematics Subject Classification:35K20, 65M06, 65M12, 82C31.
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,
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
is necessarily to be imposed for solvability of the problem. By free-energy 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 Nernst-Planck (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., ), and the membrane transport in biological ion channels (see e.g., ).
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 semi-explicit formula was derived in  for the steady-state solution; the existence and stability of the steady-state solution was established a long while ago  in the study of the steady Van Roostbroeck model in semiconductors. It has been proved by H. Gajewski and K. Gärtner in  that the solution to the drift-diffusion system converges to the thermal equilibrium state as time becomes large if the boundary conditions are in thermal equilibrium. The key-point of the proof is an energy estimate with the control of the free energy dissipation. Long time behavior was also studied in , 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 one-dimensional settings and complex three-dimensional 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 . 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 Newton-Raphson solver for the resulting linear system. In , a 3D finite element method for the system with a singular (point-like) charge outside of the concentration region was developed. A second order PNP solver was developed in  in a realistic ion-channel 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 second-order 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.,
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 long-time 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 , 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 steady-state solution. In fact, the steady-state 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  for a rigorous study of a nonlocal Poisson–Boltzmann type equation with a small dielectric parameter as derived from the two-species 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
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 Keller-Segel model (see [2, 6] and references therein). The left-hand-side in (1.6) represents the active transport of the density associated to a non-local velocity field . The potential is assumed to incorporate attractive interactions among individuals of the group, while repulsive (anti-crowding) interactions are accounted for by the diffusion in the right-hand-side.
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
The free energy contains both entropic part and the interaction part. There is a vast literature on entropic schemes for kinetic equations such as Fokker-Planck 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.
This paper is organized as follows: in Section 2, we describe the algorithm for the one dimensional case. Theoretical analysis for both semi-discrete 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
Following , we formally reformulate the system by setting
With this reformulation, an iterative algorithm may be designed as follows: given we obtain by solving the Poisson equation
We shall describe our numerical algorithm in one dimensional setting only, with one or multi-species. 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:
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
We use to approximate and to approximate , respectively. Given , we compute the potential by
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 semi-discrete approximation of the concentration satisfies
and , and .
Discretize uniformly: , so and satisfy
For the system case with the vector unknown , the PNP system in a simple form is
where is the charge of . The one-dimensional algorithm needs to be modified: (2.15) is replaced by
where , in order to solve the Poisson equation (2.18b), and each is obtained from solving
where in the formula for the potential is replaced by . For the multi-dimensional case, the algorithm can be applied in a dimension by dimension manner. We illustrate these cases by numerical examples later.
The above algorithm can be easily modified to solve the PNP equations with more physical parameters:
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.
It is also possible to generalize the schemes presented in this work to second order finite difference schemes on non-uniform 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.
The discrete concentration remains positive in time: if , then
provided the condition where
The semi-discrete free energy
A direct calculation using gives
We proceed with
since and for any and . Next we use the boundary conditions and to obtain
Putting all together leads to (2.25).
For the positivity property, we need to bound , which in multi-species 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 one-specie case . For the free energy dissipation property, we consider the semi-discrete functional of the form
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
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
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 steady-state solutions, where we notice that solutions with all initials converge to the same steady-state 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|
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 steady-state solutions. We also observe that in all cases the numerical solutions converge to the same steady-state solution for the above fixed , and , and the free free energy converges to .
3.2. Two dimensional numerical tests
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 steady-state 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.
||error in||order||error in||order|
4. Concluding remarks
In this paper, we have investigated the Poisson-Nernst-Planck 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 non-trivial flux for the potential on the boundary. We constructed simple, easy-to-implement 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 semi-discrete level. The analysis is for one and multiple ionic species and for one-dimensional problems. But numerical tests are given for both one-dimensional and two-dimensional systems, as well as for multiple species.
Although this work makes good progress in constructing and analyzing an accurate method for solving the Poisson-Nernst-Planck 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 non-zero 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.
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 13-12636.
- The same h is used for both x and y-directions.
- A. Arnold, P. Markowich and G. Toscani. On large time asymptotics for drift-diffusion-Poisson systems. Transport Theory and Statistical Physics, 29(3-5):571-58, 2000.
- J. Bedrossian, N. Rodríguez and A. L. Bertozzi. Local and global well-posedness for aggregation equations and Patlak-Keller-Segel 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 Nernst-Planck and Debye-Hückel drift-diffusion 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 Patlak-Keller-Segel 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. Three-dimensional Poisson-Nernst-Planck 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. Poisson-Nernst-Planck 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 Poisson-Boltzmann 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 Poisson-Nernst-Planck equations. J. Comput. Electron., 15:1–15, 2013
- S. Furini, F. Zerbetto and S. Cavalcanti. Application of the Poisson-Nernst-Planck theory with space-dependent 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 Poisson-Boltzmann 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 Poisson-Nernst-Planck 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 three-dimensional Poisson-Nernst-Planck 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 Poisson-Boltzmann type equations: one-dimensional 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 Poisson-Boltzmann 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 Poisson-Nernst-Planck 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 Fokker-Planck 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 Fokker-Planck 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 three-dimensional 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. Poisson-Nernst-Planck equations for simulating biomolecular diffusion-reaction 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 Poisson-Nernst-Plank 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 Poisson-Boltzmann equation. J. Phys. Chem., 94:7684–7692, 1990.
- T. Sokalski and A. Lewenstam. Application of Nernst-Planck and Poisson equations for interpretation of liquid-junction and membrane potentials in real-time and space domains. Electrochemistry Communications, 3(3):107–112, 2001.
- T. Sokalski, P. Lingenfelter and A. Lewenstam. Numerical solution of the coupled Nernst-Planck 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. Second-order Poisson-Nernst-Planck 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.