Electro-Neutral Models for Dynamic Poisson-Nernst-Planck System

Electro-Neutral Models for Dynamic Poisson-Nernst-Planck System

Zilong Song, Xiulei Cao, Huaxiong Huang Department of Mathematics and Statistics, York University and Fields Institute for Research in Mathematical Sciences, Toronto, Ontario, Canada
July 12, 2019

The Poisson-Nernst-Planck (PNP) system is a standard model for describing ion transport. In many applications, e.g., ions in biological tissues, the presence of thin boundary layers poses both modelling and computational challenges. In this paper, we derive simplified electro-neutral (EN) models where the thin boundary layers are replaced by effective boundary conditions. There are two major advantages of EN models. First of all, it is much cheaper to solve them numerically. Secondly, EN models are easier to deal with compared with the original PNP, therefore it is also easier to derive macroscopic models for cellular structures using EN models.

Even though the approach is applicable to higher dimensional cases, this paper mainly focuses on the one-dimensional system, including the general multi-ion case. Using systematic asymptotic analysis, we derive a variety of effective boundary conditions directly applicable to EN system for the bulk region. This EN system can be solved directly and efficiently without computing the solution in the boundary layer. The derivation is based on matched asymptotics, and the key idea is to bring back higher order contributions into effective boundary conditions. For Dirichlet boundary conditions, the higher order terms can be neglected and classical results (continuity of electrochemical potential) are recovered. For flux boundary conditions, however, neglecting higher contribution leads to physically incorrect solutions since they account for accumulation of ions in boundary layer. The validity of our EN model is verified by several examples and numerical computation. In particular, our EN model is much more efficient than the original PNP model when applied to the computation of membrane potential. Implemented with the Hodgkin-Huxley model, the computational time for solving the EN model is significantly reduced without sacrificing the accuracy of the solution due to the fact that it allows for relatively large mesh and time step sizes.


I Introduction

The Poisson-Nernst-Planck (PNP) system describes the transport of ions under the influence of both an ionic concentration gradient and an electric field. It is essentially a system coupling diffusion and electrostatics, and the nonlinearity comes from the drift effect of electric field on ions. Such a system and its variants have found extensive and successful applications in biological systems, in particular in the description of ion transport through cells and ion channels Gillespie and Eisenberg (2001); Horng et al. (2012). It has also been applied to many industrial fields, such as the semiconductor devices Markowich (2013) and the detection of poisonous lead by ion-selective electrode Jasielec et al. (2013).

One intriguing feature of this system is the presence of boundary layer (BL) near the boundary of concerned domain, often called Debye layer in literature. A large number of works have been devoted to the BL analysis of PNP systems. For example, singular perturbation analysis of PNP system has been carried out for narrow ion channels with certain geometric structure Singer et al. (2008); Singer and Norbury (2009). Geometric singular perturbation approach has been developed to investigate the existence and uniqueness of solutions in stationary PNP system Liu (2009); Liu and Xu (2015) as well as the effects of permanent charge and ion size Lin et al. (2013); Eisenberg et al. (2015). Recently, Wang et al. Wang et al. (2014) have tackled the steady state PNP system with arbitrary number of ion species and arbitrary valences, and have successfully reduced the asymptotic solutions to a single scalar transcendental equation.

In general, the solution consists of two parts, the BL solution in a small neighbourhood of boundary and the bulk solution in the interior region of the domain. In one-dimensional (1D) case, the leading order solution in BL can be constructed either explicitly or in integral form. Based on the BL analysis, effective continuity conditions have been proposed to connect the bulk solution and BL solution, e.g., the continuity of electro-chemical potential in Rubinstein (1990). These effective conditions have been applied to the study of steady states of 1D systems, showing the existence of multiple steady states with piecewise constant fixed charge Rubinstein (1987). One objective of this paper is to generalize the effective conditions for other boundary conditions. These conditions replace the BL region and have potential applications for deriving macroscopic models Huang et al. (2011) of bulk region in complicated structures. For example, some macroscopic continuum equations are derived in bulk region for the lens circulation Mathias (1985); Vaghefi et al. (2013), by taking into account the fluxes through membranes with an ad hoc model for the BL effect, so the fluxes calculated there might not be accurate.

The other objective of this paper is related to numerical computation of PNP systems. In addition to the BL analysis, many (conservative) numerical schemes have been developed for PNP systems, such as finite element method Gao and He (2017), finite difference scheme Flavell et al. (2014) and finite volume method Chainais-Hillairet et al. (2003); Chainais-Hillairet and Peng (2003), in one or high-dimensional spaces Lu et al. (2010); Mirzadeh and Gibou (2014). Due to the presence of BL, computation of the PNP system needs to accurately capture the behaviour of solution in BL. Since the solutions change rapidly in BL, more mesh points are needed in BL than in the bulk region to attain certain accuracy, requiring advanced techniques such as adaptive refined mesh and moving mesh Budd et al. (1996); Tang and Tang (2003). In general, computational cost is higher and development of numerical method is more demanding, especially when there are many BLs in a complicated structure. Having effective formulas/conditions to replace the BL will significantly reduce the computational time as well as the effort for developing sophisticated numerical methods, since under such framework the solutions in the bulk region can be obtained directly.

In this work, we will focus on the 1D dynamic PNP system and derive an electro-neutral (EN) system for bulk region with effective boundary conditions for several boundary conditions. In Section II, we first present the formulation for the two-ion species case and related EN models, with Dirichlet or flux boundary conditions for ion concentration and Dirichlet or Robin boundary conditions for electric potential. A more general multi-ion species model is presented afterwards. In Section III, these effective boundary conditions are validated by one steady state and two dynamic examples. In Section IV, we combine the PNP system with the Hodgkin-Huxley model and derive an EN model for neuronal axon, capturing the phenomenon of action potential efficiently. Finally conclusions and discussion of future directions are given in Section V.

Ii The electro-neutral theories

In this section, we present the electro-neutral (EN) systems with various effective boundary conditions. To introduce the main ideas, we first present the simplest PNP system, for ion specices, where the solutions and effective boundary conditions are explicit. It is followed by the general multi-ion species case.

We consider the 1D dynamic PNP system in the normalized interval ,


where and are the concentrations of the two ions with valences , and are the associated two fluxes of positive and negative ions, denotes the electric potential and is a small parameter (a combination of dielectric constant and other constants). The diffusion constants have been assumed to be 1 for simplicity, since it does not cause essential difference. Generalization to the multi-ion case with different diffusion constants will be mentioned later. We will consider various types of boundary conditions at two ends in the following subsections. For example, as in Subsection B, we can adopt Dirichlet condition for and two flux conditions at


We will also replace flux conditions by Dirichlet conditions of concentration (in Subsection C) and Dirichlet condition of by Robin-type condition (in Subsection D). The treatment will be similar at the other end . To complete the system, we also need the initial conditions for two ions. But the initial effect is not considered in this work, and we mainly limit ourselves to the large time behaviour of solutions (when BL is already present) or the case near equilibrium state.

We focus on the case when local electro-neutrality (LEN) condition in bulk region is satisfied, and there is no extra unbalanced charge present in the system/interval, or more precisely there is only unbalanced charge, here called near global electro-neutrality (NGEN) condition. We will illustrate later what kinds of boundary conditions fall in this case. These conditions can be justified in many biological applications, for example in the neuronal axon Hodgkin and Huxley (1990). Thus, in the bulk region, we assume all the functions concerned and their derivatives are , i.e.,


Then, we obtain approximately the electro-neutral condition from the first equations in (1) and more precisely we write


where and may depend on due to boundary conditions, in other words and can contain terms if boundary conditions have such terms. Substituting into second and third equation in (1) gives the EN equations


By addition and subtraction, we can also write them as


To complete this system, two effective boundary conditions are needed instead of the original three. Based on the behavior of BL solutions, we aim to derive effective conditions that connect real boundary values of (or boundary fluxes) and limit boundary values of bulk solution (or bulk fluxes). Finally we get a EN system for in the bulk region, which can be solved directly. In the following, we will always take for example, and briefly state the results for the other end.

ii.1 The leading order solution in BL

From some steady state analysis, e.g. with finite fluxes or Dirichlet conditions in Rubinstein (1990) and for Poisson-Boltzmann type equations in Lee et al. (2010) in the absence of extra unbalanced charge, there is boundary layer with thickness . Also, some numerical evidence shows that, for finite fluxes, as long as the NGEN condition is satisfied, the system has BLs near end points with all being . In this subsection, we present the leading order solutions for the PNP system. Although the solutions are well-known in literature, we give a brief derivation to be self-contained and to be more clear about the remainders.

In the BL near , we assume


and thus we set


where the argument is omitted in above functions. Then the system of equations in BL is


Integrating the second and third equations once, we get


where are the finite fluxes at . We denote as the limit values of bulk solutions at , and they should match to leading order, implying


Substituting into the first equation of (9), we get the Poisson-Boltzmann equation as leading order equation for


This can be integrated out by using and (suppose it is known here or can be expressed by known boundary conditions). Finally, we obtain


And the solutions for become


Note that in general are functions of . The composite solutions are given by


which are uniformly valid for some finite interval , say , with remainder . One can also add the contribution of BL solution near (with transform and quantities being replaced by ) to make the composite solution valid for the whole interval . Since in the bulk we have by (4), it is reasonable to expect in some intermediate region with , particularly we may choose .

ii.2 Flux boundary condition

In this subsection, we consider the case with the flux boundary conditions for two concentrations and Dirichlet condition for electric potential. More precisely, at , we have


where are given. The objective is to propose two effective boundary conditions for at based on these three functions.

To this end, we define for the EN system two fluxes


and the limit values at are denoted by respectively. Based on assumptions (3,4) in the bulk region, the two fluxes are almost the same as the two fluxes of original PNP system


where is some generic point in bulk region, say . Next, we intend to find the connection between and , or similarly between and at boundary. For this purpose, by integration of transport equations and , we immediately get


Combining these two and utilizing the composite solution (15), we obtain


where we have used the assumption that for , and by setting upper limit of integral as only exponentially small terms are neglected. In the above, may depend on . Similarly for the other flux, we obtain the relation


Physically, the quantity in above formulas is often referred to as the zeta potential Kirby (2010). To see clearly the two conditions, we carry out a linearization regarding small . In this case, they reduce to


Thus, by comparing these conditions, the total flux will almost have no difference while electric current changes, i.e.,


Physically, this means some cations/anions accumulate in the BL, and the second formula is similar to that of a capacitor. The treatment for the other end is similar, and we summarise the results below.

Proposition 1. Suppose the LEN and NGEN conditions are satisfied, and let be the given electric potential and ion fluxes at as in (16) and let be given at for original system (1), then we have the effective boundary conditions for the EN system (6)


where are defined by (17) and subscripts 0 and 1 denote quantities at and respectively.

Remark 1. Keeping the terms in the formula (24) is necessary for two reasons: first, in bulk equations (5) we have assumed an remainder so it is reasonable and consistent to bring back the terms on boundary conditions; second, neglecting the terms is physically incorrect for EN system as the solution would not be unique (e.g., can differ by a constant).
Remark 2. In this case, the fluxes can be either or , as long as the NGEN is satisfied. This means when fluxes are , we should require the fluxes are almost balanced or its integral over time satisfies


Otherwise, the solution in BL will not be anymore. For a steady state (Poisson-Boltzmann type equation in Lee et al. (2010)) with extra unbalanced charge, the solution in BL is shown to have a span of .

ii.3 Dirichlet boundary condition revisited

In this subsection, we will consider the case with Dirichlet boundary conditions for two ion concentrations. We also take left end for example, and have


The leading order effective boundary conditions for this case are well-known. With the same assumptions as previous subsection, we arrive at the same BL system, and easily get


By integration and using the matching condition, we obtain


These connection conditions are referred to as continuity of electro-chemical potential, widely adopted in literature Rubinstein (1990). And equivalently, the explicit effective boundary conditions for EN system are


As in the bulk assumption, we would like to keep the effect/terms, thus a natural question is how to bring back such effect in the effective boundary conditions for the reduced EN system. One may want to seek a general expansion to in BL and assume


The leading order solutions can be immediately written down, which are the same as those in (13,14) with replacement given by (29). However, getting the explicit expression for seems difficult. Therefore, instead we try to avoid such a solving process and intend to find the higher-order contributions directly based on leading order solution.

Now, we take for illustration, where the argument is omitted here and in the following. The second equation in BL system implies


where is some unknown flux constant. Dividing both sides by , we get


From previous subsection, we know that . Therefore, we obtain


By matching Bush (1992), let or with , we can expect that


Substituting into previous relation (33), we get from the left-hand side


and from the integral on the right-hand side


Since by definition (17), the terms automatically cancel each other (which partially verifies the correctness of matching), and we are left with


Compared with previous leading order condition (28), there is an correction term in above formula, so it can be considered as a generalization of continuity of electro-chemical potential. Treatments for the other condition and two conditions at are similar, and we summarize the results as follows.

Proposition 2.. Suppose the LEN and NGEN conditions are satisfied, and let be the given electric potential and ion concentrations at as in (26) and let be given at for original system (1), then we have the effective boundary conditions for the EN system (6)


where are defined by (17) and subscripts 0 and 1 denote quantities at and respectively.

Remark 3. We can alternatively derive an asymptotically equivalent expressions


which are Robin-type boundary conditions. Compared with (29), there are corrections in above conditions. For a special case, say at , the correction terms will be of higher order, then continuity of electro-chemical potential (29) holds with remainder .

Remark 4. In some cases, is not explicitly given, but is related to flux , so proper modification is needed. For example in biological applications, there is certain relation between flux and ion concentration across cell membrane or ion channel, such as Hodgkin-Huxley model Hodgkin and Huxley (1990) or GHK flux model Hille et al. (2001). And, in electrolyte there is the Chang-Jaffle boundary condition Chang and Jaffé (1952); Lelidis et al. (2015) or modified Chang-Jaffle condition Jasielec et al. (2013). Suppose the boundary condition is in the form , where is some given function, then we need to supplement the two conditions at with


If mixed conditions are given (e.g., one Dirichlet and one flux), then we need to combine the two relevant boundary conditions from the two propositions, e.g., if and are given, we should use and .

ii.4 Robin-type boundary condition for

In this subsection, we consider the case when Dirichlet condition of electric potential is replaced by Robin-type boundary condition. The Robin-type condition is often used to model the property of membrane or the stern layer near boundary. The previous effective conditions need to be modified since the quantity in those formulas is unknown.

ii.4.1 Dirichlet conditions for two ion concentrations

Suppose we have the boundary conditions at


where is a parameter which is assumed to be , and is some given function.

With previous assumptions, we still have (11) and (12). Integrating once and using , we obtain


or equivalently


On account of the identity , the above condition at becomes


Combining with the Robin-type condition leads to


We conclude that we have the same effective conditions as before


except that is given by (45) in this case.

Note that if , we can omit the terms in above conditions since they are not exact. If , the condition will become close to that in the Dirichlet case for . In particular, for , the term can be neglected in (45), which essentially reduces to the Dirichlet case (see (41)). If tends to infinity (not considered here), the previous BL assumptions might not be true unless , and this is left for future study.

The treatment for the right end is similar and we summarize the results below.

Proposition 3. Suppose for original system (1), the assumptions and conditions are the same as Proposition 2 except that the conditions for are replaced by


where , then we have the same effective boundary conditions for the EN system (6) as Proposition 2 except that in (38) are calculated by


where .

ii.4.2 Flux conditions

For this case, the boundary conditions at are of the form


where . The manipulation follows similar lines as before, and we summarize the results below.

Proposition 4. Suppose for original system (1), the assumptions and conditions are the same as Proposition 1 except that the conditions for are replaced by


where , then we have the same effective boundary conditions for the EN system (6) as Proposition 1 except that in (38) are determined by the nonlinear algebraic equation


where .

Remark 5. For the steady case, the above algebraic equation is the same as formula (1.23) in Lee et al. (2010), with substitution . With , the effective boundary conditions reduce to those for Dirichlet case, with .
Remark 6. The case is not considered above, since the boundary layer structure would be different. For the NGEN case considered here, it is expected that there is still a BL with thickness , and the above relation (51) implies that in BL. We postulate that in BL near ,


and it is proper to adopt the scaling and the transform instead. This seems true for Poisson-Boltzmann type equations in Lee et al. (2010) with their electro-neutral case, where boundary layer with thickness gradually disappears when becomes larger. Then, the above conditions (51,24) are still valid to leading order with new remainder , which will be verified by numerical examples in later sections.

ii.5 The general multi-ion species case

In this subsection we consider the general case with species of ions. The concentrations of ions are denoted by with valences (), where the valences are not necessarily different. The original PNP system for () and is given by


where , and are some dimensionless diffusion constants. The reduced EN system for bulk region is


where . By the LEN condition , the last concentration can be expressed by previous ones. Finally, the EN system for unknowns can be written as


where and whenever appears we should replace it by .

First, at , we consider the boundary conditions of the type


Theorem 1. Suppose the LEN and NGEN conditions are satisfied, and let be the given electric potential and ion fluxes at as in (56) and let be given at for original system (53), then for the EN system (55) we have the effective boundary conditions