An Extension of Goldman-Hodgkin-Katz Equations
by Charges from Ionic Solution and Ion Channel Protein

The Goldman-Hodgkin-Katz (GHK) equations include three equations, called the GHK flux, electric current, and voltage equations, for the calculation of ionic fluxes, electrical currents, and membrane potentials across a cell’s membrane, respectively [7, 13, 17]. They were derived by David E. Goldman and Nobel laureates Alan Lloyd Hodgkin and Bernard Katz in the 1940s [12, 14]. Since then, many studies have been done on the equations, along with a number of modified GHK equations, making them to become one of the most successful achievements in electrophysiology [4, 6, 9, 10, 15, 16, 21, 22]. As one basic tool for cell membrane study, the GHK equations have been included in most textbooks on electrophysiology (e.g. [7, 13, 17]).

Charges from an ionic solution and an ion channel protein can induce a strongly nonlinear electric filed and in turn significantly affect ionic fluxes, electrical currents, and membrane potentials. However, the GHK equations hold only for a constant electric filed since they do not consider any charge effect within an ion channel pore. Thus, they may have a low degree of approximation accuracy in the prediction of ionic fluxes, electrical currents, and voltages across a cell membrane. Even so, they are still widely applied to various ion channel studies, simulations, and model developments because they establish the connections of ionic fluxes, currents, and voltages with intracellular and extracellular concentrations and a membrane potential in algebraic expressions.

For example, in a mitochondrial dynamical system reported in [1, 5], the GHK flux equation is used to connect the concentrations from the inside of an inner mitochondrial membrane (IMM) with those from the outside of the IMM. To this end, intracellular and extracellular concentration can be treated as the variable functions of the dynamical system. Clearly, this dynamical system can be improved if the GHK flux equation is replaced by a modified GHK flux equation workable for a nonlinear electric field. Unfortunately, no such modified GHK equations existed currently.

On the other hand, lots of work was done on Poisson-Nernst-Planck ion channel (PNPic) models (e.g. see [8, 11, 18, 20] for one-dimensional (1D) PNPic models and [3, 24, 23] for recent 3D PNPic models). Although they have a higher degree of accuracy than the GHK equations in the numerical calculation of ionic fluxes, electric currents, and voltages, these PNPic models cannot be used to substitute any GHK equation since they do not provide any algebraic formula similar to a GHK equation. Hence, extending GHK equations from a constant electric field to a nonlinear electric field remains an open research problem. We intend to solve it in this work.

To do so, we revisit the derivation of GHK equations based on a typical 1D PNPic model (see (1) for its definition). We observe that one key step in the derivation of the GHK flux equation is to treat the electric field as a known coefficient function. In this way, a Nernst-Planck boundary value problem can be generated from the PNPic model and then solved analytically to obtain an algebraic formula for estimating ionic flux in terms of ionic concentration boundary values. However, expressing its solution, an ionic concentration function, in an analytical expression is difficult unless the electric field is set as a constant, as is done in the derivation of the GHK flux equation. In this work, we overcome the difficulty through introducing a new parameter, called the extension parameter, to account for the charge effects from an ionic solution and an ion channel protein. To this end, extending the GHK equations to a nonlinear electric field becomes possible.

Remarkably, we find that the extension parameter can be expressed as a definite integral of a composite function, $e^{Z_{i}u}$, over an ion channel pore interval (see (21) for details). Here, $Z_{i}$ denotes the charge number of ionic species $i$ and $u$ is a dimensionless potential function of the electric field. For a constant electric field, $u$ is simply a linear function. Thus, the extension parameter can be determined exactly. We can then reduce our extended GHK equations to the classic GHK equations. In this sense, the classic GHK equations can be regarded as special cases of our extended GHK equations. However, in the general case, we only can calculate the extension parameter approximately using a numerical quadrature. While there exist several numerical quadratures including the composite trapezoidal rule that can work for such a calculation [2], in this work, we will develop a new numerical scheme according to the composite property of the integrand function $e^{Z_{i}u}$. That is, we approximate $u$, instead of the integrand function $e^{Z_{i}u}$, as a piecewise linear interpolation function. In this way, the integral expression of the extension parameter can be rewritten as a sum of the sub-integrals that can be evaluated exactly. Clearly, our new numerical scheme can have a higher degree of numerical accuracy than the corresponding composite trapezoidal rule since its truncation error comes from a piecewise linear approximation of $u$ only. In contrast, the truncation error of the composite trapezoidal rule comes from the integrand function $e^{Z_{i}u}$.

Since a set of electrostatic potential function values is required by our new scheme, we need a methodology for its generation. One natural way for us to do so is to use a numerical solution of a 1D PNPic model. By far, several different 1D PNPic modes were developed and solved numerically (see [8, 11, 18, 20] for examples), which can be used to generate the required electrostatic potential function values. However, we could not find their software packages from the public domain. Thus, we have to develop our own software package in order to carry out an initial study on our extended GHK equations. This turns out to be one important part of this work. To do so, we will develop a nonlinear finite element iterative scheme for solving the 1D PNPic model that we have used in the construction of our extended GHK equations based on our recent 3D PNPic work [3, 24, 23]. We then will implement it in Python as a software package based on the state-of-the-art finite element library from the FEniCS project [19]. With our package, we will do numerical tests to confirm the high numerical accuracy of our specially designed numerical quadrature, to demonstrate a fast convergence rate of our finite element iterative scheme, and to show the importance of considering charge effects in the calculation of ionic fluxes.

The remaining part of this paper is organized as follows: In Section 2, we revisit the GHK equations. In Section 3, we present our extended GHK equations. In Section 4, we present our new numerical quadrature. In Section 5, we present our 1DPNPic finite element solver. In Section 6, we report numerical results. Finally, conclusions are made in Section 7.

We set a normal direction of the membrane in the $z$-axis direction of a 3D coordinate system, following what we did in our previous work [3, 24, 23]. We then define a dimensionless electrostatic potential density function, $u$, and an ionic concentration function, $c_{i}$, of species $i$ as functions of one variable, $z$, on an interval $0\leq z\leq L$ by a 1D PNPic model as follows:

	$\displaystyle-\epsilon_{s}\frac{d^{2}u}{dz^{2}}=\frac{e_{c}^{2}}{\epsilon_{0}k_{B}T}\sum_{i=1}^{n}Z_{i}c_{i}(z)+\rho(z),$	$\displaystyle 0<z<L,$			(1a)
	$\displaystyle\frac{d}{dz}{\cal D}_{i}\left[Z_{i}c_{i}(z)\frac{du(z)}{dz}+\frac{dc_{i}(z)}{dz}\right]=0,$	$\displaystyle 0<z<L,$			(1b)
	$\displaystyle u(0)=u_{0},\quad u(L)=u_{L},\quad c_{i}(0)=c_{i,0},\quad c_{i}(L)=c_{i,L},$	$\displaystyle\qquad i=1,2,\ldots,n,$			(1c)

where $\rho$ denotes a permanent charge function; ${\cal D}_{i}$ is a diffusion constant of species $i$; $n$ is the number of species in an ionic solution; $\epsilon_{s}$ denotes the water permittivity constant; $\epsilon_{0}$ is the permittivity of the vacuum; $e_{c}$ is the elementary charge; $k_{B}$ is the Boltzmann constant; $T$ is the absolute temperature; $Z_{i}$ denotes the charge number of ionic species $i$; and $u_{0}$, $u_{L}$, $c_{i,0}$, and $c_{i,L}$ are the boundary values of $u$ and $c_{i}$, respectively. Here, (1a) is called a Poisson equation, (1b) a steady Nernst-Planck equation, and (1c) the Dirichlet boundary value conditions.

Two typical selections of the interval $[0,L]$ are illustrated in Figure 1. Here, $L$ is the length of a portion of the ion channel pore, the end number $0$ is set inside the cell, and $L$ outside the cell. In Selection 1, the interval $[0,L]$ corresponds to an ion-selectivity filter of an ion channel protein. In this case, the charge function $\rho$ can be positive for an anion selectivity filter or negative for a cation selectivity filter.

The dimensionless potential $u$ can be changed to an electrostatic potential, $\Phi$, in volts by

$$\Phi(z)=\frac{k_{B}T}{e_{c}}u(z),\quad 0<z<L.$$

(2)

Here the constant $\frac{k_{B}T}{e_{c}}$ is about 0.026 volts at $T=298.15$ Kelvin.

Clearly, the steady Nernst-Planck equation (1b) can be reformulated as

$$\frac{d{\mathbf{J}}_{i}(z)}{dz}=0,\quad 0<z<L,\quad i=1,2,\ldots,n,$$

(3)

with ${\mathbf{J}}_{i}$ denoting the flux density function of ionic species $i$ as defined by

$${\mathbf{J}}_{i}(z)=-{\cal D}_{i}\left[Z_{i}c_{i}(z)\frac{du(z)}{dz}+\frac{dc_{i}(z)}{dz}\right].$$

(4)

In this work, ${\mathbf{J}}_{i}$ is a scale function because both $u$ and $c_{i}$ are the functions of one variable $z$.

We now use the above equations to derive the GHK equations by assuming that

$$\rho(z)=0,\qquad\sum_{i=1}^{n}Z_{i}c_{i}(z)=0,\quad 0<z<L.$$

In other words, none of charges are considered within an ion channel pore interval. Under this assumption, the Poisson equation (1a) becomes independent of the steady Nernst-Planck equation (1b). Thus, a two-point boundary value problem of $u$ is produced from model (1) as follows:

$$\left\{\begin{array}[]{ll}\frac{d^{2}u(z)}{dz^{2}}=0,&0<z<L\\ u(0)=u_{0},\quad u(L)=u_{L}.&\end{array}\right.$$

The solution $u$ of the above boundary value problem can be easily found in the expression

$$u(z)=\frac{u_{L}-u_{0}}{L}z+u_{0},\quad 0\leq z\leq L,$$

(5)

resulting in a constant electric field in the form

$$\frac{du(z)}{dz}=-\frac{V_{m}}{L},$$

(6)

where $V_{m}$ is defined by

$$V_{m}=u_{0}-u_{L},$$

(7)

which will be referred to as a membrane potential (or a voltage across a cell membrane).

From (3) it implies that the flux function ${\mathbf{J}}_{i}(z)$ is either 0 or a nonzero constant. In the case in which ${\mathbf{J}}_{i}=0$, we use (4) to get

$$Z_{i}c_{i}(z)\frac{du(z)}{dz}+\frac{dc_{i}(z)}{dz}=0,\quad 0<z<L.$$

By separating the functions $u$ and $c_{i}$, the above equation is rewritten as

$$\frac{1}{c_{i}}\frac{dc_{i}}{dz}=-Z_{i}\frac{du}{dz}.$$

Integrating over the interval $(0,L)$ on the both sides of the above equation, we obtain the Nernst equation:

$$V_{m}=\frac{1}{Z_{i}}\ln\frac{c_{i}(L)}{c_{i}(0)},\quad i=1,2,\ldots,n.$$

(8)

We next consider the case in which ${\mathbf{J}}_{i}$ is a constant. Since the electric field $\frac{du}{dz}$ is a constant, as given in (6), we can use (4) to construct a two-point boundary value problem for determining the concentration $c_{i}$ of ionic species $i$ as follows:

$$\displaystyle\left\{\begin{array}[]{ll}\displaystyle\frac{dc_{i}(z)}{dz}-Z_{i}\frac{V_{m}}{L}c_{i}(z)=-\frac{{\mathbf{J}}_{i}}{{\cal D}_{i}},&0<z<L,\\ c_{i}(0)=c_{i,0},\quad c_{i}(L)=c_{i,L}.&\end{array}\right.$$

(9)

Since $\frac{{\mathbf{J}}_{i}}{{\cal D}_{i}}$ is a constant, by the integrating factor method, $c_{i}$ can be found in the expression

$$c_{i}(z)=\frac{1}{\theta(z)}\left[-\frac{{\mathbf{J}}_{i}}{{\cal D}_{i}}\int_{Z1}^{z}\theta(z)dz+c_{i,0}\right],$$

(10)

where $\theta(z)$ is an exponential integrating factor, which can be found as follows:

$$\theta(z)=e^{-\int_{0}^{z}Z_{i}\frac{V_{m}}{L}dz}=e^{-Z_{i}\frac{V_{m}}{L}z},$$

where $V_{m}$ has been defined in (7). Applying the above expression to (10), we get

$$c_{i}(z)=e^{Z_{i}\frac{V_{m}}{L}z}\left[\frac{{\mathbf{J}}_{i}L}{{\cal D}_{i}Z_{i}V_{m}}\left(e^{-Z_{i}\frac{V_{m}}{L}z}-1\right)+c_{i,0}\right].$$

(11)

We then set $z=L$ on the both sides of (11) to get a linear equation of ${\mathbf{J}}_{i}$:

$$c_{i}(L)=e^{Z_{i}V_{m}}\left[\frac{{\mathbf{J}}_{i}L}{{\cal D}_{i}Z_{i}V_{m}}\left(e^{-Z_{i}V_{m}}-1\right)+c_{i,0}\right].$$

Solving the above equation for ${\mathbf{J}}_{i}$, we obtain an expression of ${\mathbf{J}}_{i}$ as follows:

$${\mathbf{J}}_{i}=Z_{i}{\cal D}_{i}\left(\frac{V_{m}}{L}\right)\frac{c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}}{1-e^{-Z_{i}V_{m}}},\qquad i=1,2,\ldots,n,$$

(12)

which can also be rewritten as

$${\mathbf{J}}_{i}=Z_{i}{\cal D}_{i}\left(\frac{V_{m}}{L}\right)\frac{c_{i,L}-c_{i,0}e^{Z_{i}V_{m}}}{1-e^{Z_{i}V_{m}}}.$$

(13)

The expression (12) or (13) is referred to as the GHK flux equation in the literature.

The electrical current ${\cal I}_{i}$ of species $i$ can be derived from the multiplication of ${\mathbf{J}}_{i}$ by the charge $Z_{i}e_{c}$ of species $i$:

$${\cal I}_{i}=Z_{i}e_{c}{\mathbf{J}}_{i},\quad i=1,2,\ldots,n.$$

Using the GHK flux equation (12), we immediately derive the GHK electric current equation:

$${\cal I}_{i}=Z_{i}^{2}e_{c}{\cal D}_{i}\left(\frac{V_{m}}{L}\right)\frac{c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}}{1-e^{-Z_{i}V_{m}}}.$$

(14)

Adding the above currents together, we obtain a net electrical current, denoted by ${\cal I}$, as follows:

$${\cal I}=\sum_{i=1}^{n}{\cal I}_{i}=e_{c}\frac{V_{m}}{L}\sum_{i=1}^{n}Z_{i}^{2}{\cal D}_{i}\frac{c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}}{1-e^{-Z_{i}V_{m}}}.$$

(15)

If the net current is zero, setting ${\cal I}=0$, we derive the GHK voltage equation:

$$\sum_{i=1}^{n}Z_{i}^{2}{\cal D}_{i}\frac{c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}}{1-e^{-Z_{i}V_{m}}}=0.$$

(16)

A solution $V_{m}$ of the above equation is called the GHK potential.

For example, for a solution with sodium (Na⁺ with $Z_{1}=1$), potassium (K⁺ with $Z_{2}=1$), and chloride (Cl^- with $Z_{3}=-1$) ions ($n=3$) [7, Eq. (2.3)] and [17, Eq. (2.7.2)], we can solve the equation (16) for $V_{m}$ to get the GHK potential in the expression:

$$V_{m}=-\ln\frac{{\cal D}_{1}c_{1,0}+{\cal D}_{2}c_{2,0}+{\cal D}_{3}c_{3,L}}{{\cal D}_{1}c_{1,L}+{\cal D}_{2}c_{2,L}+{\cal D}_{3}c_{3,0}},$$

where $c_{1,0},c_{2,0}$ and $c_{3,0}$ are the concentration values of sodium, potassium, and chloride ions at $z=0$ while $c_{1,L},c_{2,L}$ and $c_{3,L}$ at $z=L$, respectively.

A membrane potential, $V$, in volts is defined by

$$V=\Phi(0)-\Phi(L).$$

Using (2) and (7), we can derive $V$ from $V_{m}$ by

$$V=\frac{k_{B}T}{e_{c}}V_{m}.$$

(17)

Using (12) and (17), we can obtain another expression of the GHK flux equation as follows:

$${\mathbf{J}}_{i}=Z_{i}{\cal D}_{i}\left(\frac{e_{c}}{k_{B}TL}\right)V\frac{c_{i,0}-c_{i,L}e^{-\frac{Z_{i}e_{c}}{k_{B}T}V}}{1-e^{-\frac{Z_{i}e_{c}}{k_{B}T}V}},\quad i=1,2,\ldots,n.$$

(18)

As a special case, when $c_{i,0}=c_{i,L}=c_{i}^{b}$, the above expression is simplified as

$${\mathbf{J}}_{i}=Z_{i}{\cal D}_{i}c_{i}^{b}\left(\frac{e_{c}}{k_{B}TL}\right)V,$$

indicating that the flux ${\mathbf{J}}_{i}$ has a linear relationship with the membrane potential $V$. In other words, the GHK flux equation describes a nonlinear relationship of ${\mathbf{J}}_{i}$ with $V$ only if $c_{i,0}\neq c_{i,L}$.

In this section, we extend the GHK equations from a constant electric field to a nonlinear electric field based on the 1D PNPic model (1) with a nonzero permanent charge function $\rho$. Because of (3), we still have that the flux ${\mathbf{J}}_{i}$ of species $i$ is either zero or a nonzero constant. Since the case of ${\mathbf{J}}_{i}$ being zero can be treated in the same way as what is done in the previous section, we only consider the case of ${\mathbf{J}}_{i}$ being a nonzero constant in this section.

Similarly to what is done in the derivation of GHK flux equation, we treat the nonlinear electric field $du/dz$ as a known coefficient function of the steady Nernst-Planck equation (1b). In this way, we can separate (1b) from the Poisson equation (1a) and then generate a two-point boundary value problem for each ionic concentration function, $c_{i}$, from the 1D PNPic model (1) as follows:

$$\displaystyle\left\{\begin{array}[]{ll}\displaystyle\frac{dc_{i}}{dz}+Z_{i}\frac{du}{dz}c_{i}(z)=-\frac{{\mathbf{J}}_{i}}{{\cal D}_{i}},&0<z<L,\\ c_{i}(0)=c_{i,0},\quad c_{i}(L)=c_{i,L}.&\end{array}\right.$$

By the integrating factor method, the solution $c_{i}$ of the above problem can be found in the form

$$c_{i}(z)=\frac{1}{\theta(z)}\left[-\frac{{\mathbf{J}}_{i}}{{\cal D}_{i}}\int_{0}^{z}\theta(\xi)d\xi+c_{i,0}\right],\quad 0<z<L,$$

(19)

where $\theta(z)$ denotes an exponential integrating factor, which can be found as shown below:

$$\displaystyle\theta(z)=e^{\int_{0}^{z}Z_{i}\frac{du(\xi)}{d\xi}d\xi}=e^{Z_{i}[u(z)-u_{0}]}.$$

Applying the above expression to (19), we obtain an integral expression of $c_{i}$ as follows:

$$c_{i}(z)=e^{-Z_{i}[u(z)-u_{0}]}\left[-\frac{{\mathbf{J}}_{i}}{{\cal D}_{i}}\int_{0}^{z}e^{Z_{i}[u(\xi)-u_{0}]}d\xi+c_{i,0}\right].$$

Setting $z=L$ on the both sides of the above expression, we get a linear equation of ${\mathbf{J}}_{i}$:

$$c_{i,L}=e^{-Z_{i}[u_{L}-u_{0}]}\left[-\frac{{\mathbf{J}}_{i}}{{\cal D}_{i}}\int_{0}^{L}e^{Z_{i}[u(\xi)-u_{0}]}d\xi+c_{i,0}\right].$$

Solving the above equation for ${\mathbf{J}}_{i}$, we get an extension of the GHK flux equation in the form

$${\mathbf{J}}_{i}={\cal D}_{i}\frac{c_{i,0}-c_{i,L}e^{Z_{i}[u_{L}-u_{0}]}}{\int_{0}^{L}e^{Z_{i}[u(z)-u_{0}]}dz},\quad i=1,2,\ldots,n,$$

(20)

Clearly, the integral $\int_{0}^{L}e^{Z_{i}[u(z)-u_{0}]}dz$ can be reformulated as

$$\int_{0}^{L}e^{Z_{i}[u(z)-u_{0}]}dz=e^{-Z_{i}u_{0}}\int_{0}^{L}e^{Z_{i}u(z)}dz.$$

Note that the integral $\int_{0}^{L}e^{Z_{i}u(z)}dz$ is a constant that collects all the values of potential function $u$ over the ion channel interval $[0,L]$. Hence, we can use it to define an extension parameter, $\alpha_{i}$, by

$$\alpha_{i}=\int_{0}^{L}e^{Z_{i}u(z)}dz,\quad i=1,2,\ldots,n.$$

(21)

For clarity, we denote the extended GHK flux ${\mathbf{J}}_{i}$ of (20) by the new notation ${\mathbf{J}}_{i}^{E}$ and use the extension parameter $\alpha_{i}$ to reformulate it as follows:

$${\mathbf{J}}_{i}^{E}={\cal D}_{i}\frac{e^{Z_{i}u_{0}}}{\alpha_{i}}\left[c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}\right],\quad i=1,2,\ldots,n,$$

(22)

where $V_{m}$ and $\alpha_{i}$ have been defined in (7) and (21), respectively.

Following what is done in (14), we can obtain an extened GHK current equation as follows:

$${\cal I}_{i}^{E}=e_{c}Z_{i}{\cal D}_{i}\frac{e^{Z_{i}u_{0}}}{\alpha_{i}}\left[c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}\right],\quad i=1,2,\ldots,n,$$

(23)

where ${\cal I}_{i}^{E}$ denotes an electrical current of species $i$.

Adding all ${\cal I}_{i}^{E}$ together, we obtain a net electrical current, ${\cal I}^{E}$, in the expression:

$$\quad{\cal I}^{E}=e_{c}\sum_{i=1}^{n}Z_{i}{\cal D}_{i}\frac{e^{Z_{i}u_{0}}}{\alpha_{i}}\left[c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}\right].$$

(24)

Similar to what is done in (16), we can set ${\cal I}^{E}=0$ to obtain an extended GHK voltage equation as follows:

$$\sum_{i=1}^{n}Z_{i}{\cal D}_{i}\frac{e^{Z_{i}u_{0}}}{\alpha_{i}}\left[c_{i,0}-c_{i,L}e^{-Z_{i}V_{m}}\right]=0.$$

(25)

A solution of the above equation gives an extended ned GHK potential.

As an example, from the equation of (25) we can get an extended ned GHK potential, denoted by $V_{m}^{E}$, for an ionic solution with sodium (Na⁺ with $Z_{1}=1$) and chloride (Cl^- with $Z_{2}=-1$) ions ($n=2$) in the expression

$$V_{m}^{E}=-\ln\frac{\varpi+\sqrt{\varpi^{2}+4\frac{\alpha_{1}}{\alpha_{2}}{\cal D}_{1}{\cal D}_{2}e^{-2u_{0}}c_{1,L}c_{2,L}}}{2{\cal D}_{1}c_{1,L}},$$

(26)

where $c_{1,L}$ and $c_{2,L}$ denote the concentration values of sodium and chloride ions at $z=L$, respectively; ${\cal D}_{1}$ and ${\cal D}_{2}$ are the diffusion constants of sodium and chloride ions, respectively; $\alpha_{1}$, $\alpha_{2}$, and $\varpi$ are defined by

$$\alpha_{1}=\int_{0}^{L}e^{u(z)}dz,\quad\alpha_{2}=\int_{0}^{L}e^{-u(z)}dz,$$

and

$$\varpi={\cal D}_{1}c_{1,0}-\frac{\alpha_{1}}{\alpha_{2}}{\cal D}_{2}e^{-2u_{0}}c_{2,0}.$$

Here $\varpi$ has been assumed to be nonnegative and we have selected the positive square root in (26) to ensure that the extended ned GHK potential $V_{m}^{E}$ is well-defined.

In this section, we present a numerical quadrature scheme for computing the extension parameter $\alpha_{i}$ approximately. In this scheme, we assume that a set of interpolation points, ${\cal S}_{m}$, is given by

$${\cal S}_{m}=\{(z_{j},u_{j})\;|\;j=0,1,2,\ldots,m\},\quad m\geq 1,$$

(27)

where $z_{j}$ denotes the $j$-th partition number of the ion channel pore interval $[0,L]$ satisfying

$$0=z_{0}<z_{1}<z_{2}<\ldots<z_{m-1}<z_{m}=L,$$

(28)

and $u_{j}$ denotes a numerical value of an electrostatic potential function, $u$, at $z=z_{j}$ for $j=0,1,2,\ldots,m$. To simplify the presentation of our numerical scheme, we set $z_{j}=jh$ with the mesh size $h=L/m$.

When $u$ is nonlinear, we use the interpolation point set ${\cal S}_{m}$ to construct a piecewise linear interpolating function, $\hat{u}$, of $u$ in the expression

$$\hat{u}(z)=a_{j}z+b_{j}\quad\mbox{for }z_{j-1}\leq z\leq z_{j},$$

(29)

where $a_{j}=(u_{j}-u_{j-1})/h$ and $b_{j}=(z_{j}u_{j-1}-z_{j-1}u_{j})/h$ for $j=1,2,\ldots,m$. That is, $\hat{u}$ is a linear function within sub-interval $[z_{j-1},z_{j}]$. We then estimate the extension parameter $\alpha_{i}$ approximately as shown below:

	$\displaystyle\alpha_{i}$	$\displaystyle=$	$\displaystyle\int_{0}^{L}e^{Z_{i}u(z)}dz=\sum_{j=1}^{m}\int_{z_{j-1}}^{z_{j}}e^{Z_{i}u(z)}dz$
		$\displaystyle\approx$	$\displaystyle\sum_{j=1}^{m}\int_{z_{j-1}}^{z_{j}}e^{Z_{i}\hat{u}(z)}dz=\sum_{j=1}^{m}\int_{z_{j-1}}^{z_{j}}e^{Z_{i}(a_{j}z+b_{j})}dz$
		$\displaystyle=$	$\displaystyle\frac{1}{Z_{i}}\sum_{j=1}^{m}\frac{e^{Z_{i}b_{j}}}{a_{j}}\left[e^{Z_{i}a_{j}z_{j}}-e^{Z_{i}a_{j}z_{j-1}}\right]$
		$\displaystyle=$	$\displaystyle\frac{h}{Z_{i}}\sum_{j=1}^{m}\frac{e^{Z_{i}u_{j}}-e^{Z_{i}u_{j-1}}}{u_{j}-u_{j-1}},$

where we have used the partition formula $z_{j}=z_{j-1}+h$. For clarity, the above approximate value of $\alpha_{i}$ is denoted by the notation $\hat{\alpha}_{i}$. That is,

$$\hat{\alpha}_{i}=\frac{h}{Z_{i}}\sum_{j=1}^{m}\frac{e^{Z_{i}u_{j}}-e^{Z_{i}u_{j-1}}}{u_{j}-u_{j-1}}.$$

(30)

This gives the formula for computing the extended GHK flux ${\mathbf{J}}_{i}^{E}$ approximately as follows:

$${\mathbf{J}}_{i}^{E}\approx{\cal D}_{i}\frac{e^{Z_{i}u_{0}}}{\hat{\alpha}_{i}}\left[c_{i}(0)-c_{i}(L)e^{-Z_{i}V_{m}}\right],\quad i=1,2,\ldots,n,$$

(31)

where $V_{m}$ has been defined in (7). Specifically, when $m=1$, we have

$$h=L,\quad z_{0}=0,\quad z_{1}=L,\quad u_{0}=u(0),\quad u_{1}=u(L),$$

and

$$\hat{\alpha}_{1}=\frac{h}{Z_{i}}\frac{e^{Z_{i}u_{1}}-e^{Z_{i}u_{0}}}{u_{1}-u_{0}},$$

from which (31) can be reduced to the GHK flux equation (12). In this sense, our extended GHK flux equation contains the classic GHK flux equation as a special case (i.e. only using the two boundary values $u_{0}$ and $u_{L}$ of $u$).

In order to generate an interpolation date set of (27), we present a nonlinear finite element iterative scheme for solving the 1DPNPic model (1) in this section.

We start with a finite element approximation of the 1DPNPic model (1).

Let $V$ denote a linear finite element function space based on the mesh partition (28) of the interval $[0,L]$. Here each function of $V$ is linear within each subinterval of the partition and continuous on the interval $[0,L]$. We define a subspace, $V_{0}$, of $V$ by

$$V_{0}=\{v\in V\;|\;v(0)=0,\quad v(L)=0\}.$$

In the finite element method, $V$ and $V_{0}$ are referred to as the trial and test function spaces, respectively. Multiplying each equation of the system (1) by a test function, $v$, and integrating the second-order terms by parts, we can approximate the 1DPNPic model (1) as a system of nonlinear finite element variational problems: Find $u\in V$ and $c_{i}\in V$ satisfying the Dirichlet boundary value conditions $u(0)=u_{0}$, $u(L)=u_{L}$, $c_{i}(0)=c_{i,0}$, and $c_{i}(L)=c_{i,L}$ such that

	$\displaystyle\epsilon_{s}\int_{0}^{L}\frac{du}{dz}\frac{dv}{dz}dz-\frac{e_{c}^{2}}{\epsilon_{0}k_{B}T}\int_{0}^{L}\left(\sum_{j=1}^{n}Z_{j}c_{j}\right)vdz=\int_{0}^{L}\rho(z)vdz\quad\forall v\in V_{0},$		(32a)
	$\displaystyle\int_{0}^{L}{\cal D}_{i}\left(Z_{i}c_{i}\frac{du}{dz}+\frac{dc_{i}}{dz}\right)\frac{dv_{i}}{dz}dz=0\quad\forall v_{i}\in V_{0},\quad i=1,2,\ldots,n.$		(32b)

We next present a damped iterative scheme for solving the above nonlinear finite element system.

Let $u^{(k)}$ and $c_{i}^{(k)}$ denote the $k$th iterates of $u$ and $c_{i}$, respectively. When the initial iterates $u^{(0)}$ and $c_{i}^{(0)}$ are given, for $k=0,1,2,\ldots,$ we define the updates $u^{(k+1)}$ and $c_{i}^{(k+1)}$ by

	$\displaystyle c_{i}^{(k+1)}$	$\displaystyle=$	$\displaystyle c_{i}^{(k)}+\omega({p}_{i}-c_{i}^{(k)}),\quad i=1,2,\ldots,n,$		(33a)
	$\displaystyle u^{(k+1)}$	$\displaystyle=$	$\displaystyle u^{(k)}+\omega({q}-u^{(k)}),$		(33b)

where $\omega$ is a damping parameter between 0 and 1, ${p}_{i}$ is a solution of the linear finite element variational problem: Find $p_{i}\in V$ satisfying $p_{i}(0)=c_{i,0}$ and $p_{i}(L)=c_{i,L}$ for $i=1,2,\ldots,n$ such that

$$\int_{0}^{L}{\cal D}_{i}\left(Z_{i}p_{i}\frac{du^{(k)}}{dz}+\frac{dp_{i}}{dz}\right)\frac{dv_{i}}{dz}dz=0\quad\forall v_{i}\in V_{0},$$

(34)

and ${q}$ is a solution of the linear finite element variational problem: Find $q\in V$ satisfying $q(0)=u_{0}$ and $q(L)=u_{L}$ such that

$$\epsilon_{s}\int_{0}^{L}\frac{dq}{dz}\frac{dv}{dz}dz=\frac{e_{c}^{2}}{\epsilon_{0}k_{B}T}\int_{0}^{L}\left(\sum_{j=1}^{n}Z_{j}c_{j}^{(k+1)}\right)vdz+\int_{0}^{L}\rho(z)vdz\quad\forall v\in V_{0}.$$

(35)

Note that we have separate the nonlinear system (32) into $n+1$ independent linear equations by substituting the current iterate $u^{(k)}$ to $u$ in (32a) and the updates $c_{i}^{(k+1)}$ to $c_{i}$ in (32b). By default, the initial iterates $u^{(0)}$ and $c_{i}^{(0)}$ are set by

$$u^{(0)}=\frac{u_{L}-u_{0}}{L}z+u_{0},\quad c_{i}^{(0)}=0,\quad i=1,2,\ldots,n,$$

(36)

since the above selection can lead to the GHK flux equation as shown in Section 2. We terminate the iterative process whenever the following iteration termination rule holds:

$$\|u^{(k+1)}-u^{(k)}\|<\epsilon,\quad\max_{1\leq i\leq n}\|{c}^{(k+1)}_{i}-{c}^{(k)}_{i}\|<\epsilon,$$

(37)

where $\epsilon$ is a tolerance (e.g. $\epsilon=10^{-5}$) and $\|\cdot\|$ is defined by

$$\|v\|=\sqrt{\int_{0}^{L}|v(z)|^{2}dz}\quad\mbox{for }v\in V.$$

We implemented the damped iterative scheme (33) in Python as a program package based on the state-of-the-art finite element library from the FEniCS project [19]. In numerical tests, we used a piecewise constant permanent charge function, $\rho$, in the expression

$$\rho(z)=\left\{\begin{array}[]{ll}Q_{1},&0<z\leq 10,\\ Q_{2},&10<z\leq 20,\\ Q_{3},&20<z\leq 30,\\ Q_{4},&30<z\leq 40,\\ \end{array}\right.$$

(38)

where $Q_{i}$ is a net charge within the $i$th sub-interval of the ion channel interval $[0,L]$ with $L=40$ for $i=1,2,3,4$. When all $Q_{i}$ are equal to $Q$, $\rho$ becomes a constant function. We did three numerical tests, called Tests 1, 2, and 3, using the following values of $Q$ and $Q_{i}$:

Test 1

Set $Q_{i}=Q$ for $i=1,2,3,4$ with $Q=10$.

Test 2

Set $Q_{i}=Q$ for $i=1,2,3,4$ with $Q=-20$.

Test 3

Set $Q_{1}=0,Q_{2}=20,Q_{3}=10,$ and $Q_{4}=0$.

Here Tests 1 and 2 are used to study the extended GHK flux equation and 1DPNPic model within an anion-selectivity filter and a cation-selectivity filter, respectively. In Test 3, the interval $[0,L]$ contains a selectivity filter (i.e. $10<z\leq 30$) and two solution regions (i.e. $(0,10]$ and $(30,40]$), and $\rho$ is a non-constant positive charge function within the filter and zero within the two neutral solution regions. See Figure 1(a) for an illustration of the ion channel interval $[0,L]$ used in Tests 1 and 2 and Figure 1(b) for that used in Test 3.

We also did numerical tests without considering any permanent charge effect as follows:

Test 4

Set $\rho=0$ and $c_{i,0}\neq c_{i,L}$ with $c_{i,0}=0.1$ and $c_{i,L}=0.5$ for $i=1,2$.

Test 5

Set $\rho=0$ and $c_{i,0}=c_{i,L}$ with $c_{i,0}=0.1$ and $c_{i,L}=0.1$ for $i=1,2$.

In these numerical tests, the units of length and concentration were set as angstroms and moles per liter (mol/L), respectively; a mesh partition of (28) was set with $m=256$ and $L=40$, which gave the mesh size $h\approx 0.156$; a solution of two ionic species (i.e. $n=2$) was used with $c_{1}$ denoting a concentration of cations and $c_{2}$ a concentration of anions; the boundary values of (1c) were set as $u_{0}=-2$, $u_{L}=2,$ $c_{i,0}=0.1$, and $c_{i,L}=0.5$ for $i=1,2$. We also fixed the parameters $\epsilon_{s}=80$, $Z_{1}=1$, ${\cal D}_{1}=0.133$ (for Na⁺), $Z_{2}=-1$, and ${\cal D}_{2}=0.203$ (for Cl^-). The initial iterates were set by (36). Each finite element equation of (34) and (35) were solved by the Gaussian elimination method. The numerical tests were done on our iMac computer with one 4.2 GHz Intel core i7 processor and 64 GB memory. The numerical results are reported in Table 1 and Figures 2 and 3.

Table 1 reports the flux values ${\mathbf{J}}_{1}^{E}$ and ${\mathbf{J}}_{2}^{E}$ calculated by our extended GHK flux equations (31) and their relative errors with the ionic fluxes ${\mathbf{J}}_{1}$ and ${\mathbf{J}}_{2}$ calculated by the GHK flux equation (12). The relative error values can be large up to about 7, showing the importance of considering charge effects. Table 1 also reports the convergence and performance of our nonlinear iterative scheme (33) in terms of the number of iterations and computer CPU time. Here the damping parameter $\omega$ was selected optimally (by numerical experiments) in the sense that it enabled the number of iterations satisfying the iteration termination rule (37) to become the smallest. The computer CPU time was counted from the starting of an implementation to the completion of an iteration process. From the table it can be seen that our finite element iterative scheme (33) has a fast convergence rate and runs efficiently. All the tests only took only 0.33 seconds or less in CPU time.

Figure 2 displays the electrostatic potential functions and ionic concentration functions generated by our finite element program package for Tests 1, 2, and 3. In the case of Test 1, from Figure 2(a) it can be seen that anions have dominated cations due to a positive net charge within the selectivity filter of an anion channel protein (such as a chloride channel). On the other hand, when a negative net charge is set in the selectivity filter of an ion channel protein, cations become to dominate anions as shown in Figure 2(b). This validates the selectivity of a cation channel (such as a potassium channel). In Test 3, we used a nonnegative piecewise constant function to mimic a charge distribution over a portion of an ion channel pore that contains one anion selectivity filter and two neutral solution regions, resulting in a dominating concentration of anions as shown in Figure 2(c). These test results indicate that our finite element program package is a valuable tool in the study of the ion selectivity property of an ion channel protein.

Interestingly, the numerical results for Tests 4 and 5 reported in Table 1 and Figure 3 show that the boundary values of ionic concentrations $c_{i}$ can significantly affect the electrostatic potential and ionic distributions when $\rho=0$ (i.e. ignoring the charge effects from an ion channel protein). When the boundary values are different, as was done in Test 4, our GHK flux equation still produced different ionic fluxes from the classic GHK flux equation. However, with the same boundary values, we got a linear potential function, $u$, of $z$ and two constant ionic concentrations $c_{i}$ for $i=1,2$, as shown in Figure 3(b). As a result, our extended GHK flux equation produced the same ionic fluxes as the classic GHK equation.