Geometrical effects on nonlinear electrodiffusion in cell physiology

TL;DR

This paper develops new electrical laws based on nonlinear electro-diffusion theory to understand how local geometrical features like curvature influence the electrical properties of cells, with implications for nano-pipette design and neural signal analysis.

Contribution

It introduces an asymptotic approximation for electro-diffusion in complex geometries, revealing how charge concentration depends on local curvature and shape.

Findings

Charge concentrates at cusp-shaped features as charge increases

Derived asymptotic solutions for steady-state electro-diffusion in complex geometries

Insights into voltage variations in neurons with heterogeneous structures

Abstract

We report here new electrical laws, derived from nonlinear electro-diffusion theory, about the effect of the local geometrical structure, such as curvature, on the electrical properties of a cell. We adopt the Poisson-Nernst-Planck (PNP) equations for charge concentration and electric potential as a model of electro-diffusion. In the case at hand, the entire boundary is impermeable to ions and the electric field satisfies the compatibility condition of Poisson's equation. We construct an asymptotic approximation for certain singular limits to the steady-state solution in a ball with an attached cusp-shaped funnel on its surface. As the number of charge increases, they concentrate at the end of cusp-shaped funnel. These results can be used in the design of nano-pipettes and help to understand the local voltage changes inside dendrites and axons with heterogenous local geometry.

Equations

$$\displaystyle\int\limits_{\Omega}q(\mbox{\boldmath$x$})\,d\mbox{\boldmath$x$}=N.$$

$$Q=zeN,$$

$$\displaystyle D\left[\Delta\rho(\mbox{\boldmath$x$},t)+\frac{ze}{kT}\nabla\left(\rho(\mbox{\boldmath$x$},t)\nabla\phi(\mbox{\boldmath$x$},t)\right)\right]=$$

$$\displaystyle D\left[\frac{\partial\rho(\mbox{\boldmath$x$},t)}{\partial n}+\frac{ze}{kT}\rho(\mbox{\boldmath$x$},t)\frac{\partial\phi(\mbox{\boldmath$x$},t)}{\partial n}\right]=$$

$$\displaystyle\rho(\mbox{\boldmath$x$},0)=$$

$$\displaystyle\Delta\phi(\mbox{\boldmath$x$},t)=$$

$$\displaystyle\frac{\partial\phi(\mbox{\boldmath$x$},t)}{\partial n}=$$

$$\displaystyle\sigma(\mbox{\boldmath$x$},t)=\frac{Q}{\varepsilon\varepsilon_{0}|\partial\Omega|}.$$

$$\displaystyle\rho(\mbox{\boldmath$x$})=N\frac{\exp\left\{-\displaystyle\frac{ze\phi(\mbox{\boldmath$x$})}{kT}\right\}}{\displaystyle\int_{\Omega}\exp\left\{-\displaystyle\frac{ze\phi(\mbox{\boldmath$x$})}{kT}\right\}\,d\mbox{\boldmath$x$}},$$

$$\displaystyle\Delta\phi(\mbox{\boldmath$x$})=-\frac{zeN\exp\left\{-\displaystyle\frac{ze\phi(\mbox{\boldmath$x$})}{kT}\right\}}{\varepsilon\varepsilon_{0}{\displaystyle\int_{\Omega}\exp\left\{-\displaystyle\frac{ze\phi(\mbox{\boldmath$x$})}{kT}\right\}\,d\mbox{\boldmath$x$}}}{.}$$

$$\displaystyle\frac{\partial\phi({\mbox{\boldmath$x$}})}{\partial n}=-\frac{Q}{\varepsilon\varepsilon_{0}|\partial\Omega|},$$

$$\displaystyle u(\mbox{\boldmath$x$})=\displaystyle\frac{ze\phi({\mbox{\boldmath$x$}})}{kT},\quad\lambda=\frac{(ze)^{2}N}{\varepsilon\varepsilon_{0}kT},$$

$$\displaystyle\Delta u(\mbox{\boldmath$x$})=$$

$$\displaystyle\frac{\partial u(\mbox{\boldmath$x$})}{\partial n}=-\frac{\lambda}{|\partial\Omega|}\hskip 5.0pt\mbox{for}\ \mbox{\boldmath$x$}\in\partial\Omega.$$

$$\displaystyle-\Delta\tilde{u}(\mbox{\boldmath$x$})=$$

$$\displaystyle\frac{\partial\tilde{u}(\mbox{\boldmath$x$})}{\partial n}=$$

$$w=w(z)=\frac{z-\alpha }{1-\alpha z},$$

$$\alpha =-1-\sqrt{\epsilon }+O(\epsilon ).$$

$${\Delta }_{w}v(w)=$$

$$=$$

$$\frac{\partial u(z)}{\partial n_z}={-\frac{\partial v(w)}{\partial \theta }|}_{w=-1}\frac{\partial \theta }{\partial Y},$$

$$i{e}^{i\theta }\frac{\partial \theta }{\partial Y}=w^{\prime }(z)=\frac{1-{\alpha }^2}{(1-\alpha z{)}^2}.$$

$$\frac{\partial v(w)}{\partial n}=-\frac{\lambda \sqrt{\epsilon }}{2|\partial \Omega |}\text{mbox}forw\in \partial {\Omega }_{w,a}.$$

$$v^{\prime \prime }+\frac{4\epsilon }{|e^{i\theta }-1-e^{i\theta }\sqrt{\epsilon }{|}^4}\exp \left\{-v(e^{i\theta })\right\}=$$

$$v^{\prime }(c\sqrt{\epsilon })=$$

$$v^{\prime }(\pi )=$$

$$\tilde{\Delta }u=u(N)-u(C)=v(c\sqrt{\epsilon })-v(\pi ).$$

$$y_{\text{mbox}outer}(\theta )=M\theta +M^{\prime },$$

$$\frac{{\epsilon }^2}{|e^{i\theta }-1-e^{i\theta }\sqrt{\epsilon }{|}^4}=\frac{1}{(1+{\xi }^2{)}^2}+O(\sqrt{\epsilon }).$$

$$Y^{\prime \prime }(\xi )+\frac{4}{(1+{\xi }^2{)}^2}\exp \left\{-Y(\xi )\right\}=0,$$

$$Y(\xi )=A+B_2{\xi }^2+B_4{\xi }^4+\dots ,$$

$$B_2=-2e^{-A}.$$

$$Y^{\prime \prime }(\xi )+\frac{4e^{-A}}{(1+{\xi }^2{)}^2}=0$$

$$Y_{appr}(\xi )=-2\xi e^{-A}\arctan \xi .$$

Full text

1 Introduction

Electro-diffusion is the process by which the motion of ions in solution is driven by two physical forces: thermal motion, which is diffusion, and the electric field. The difficulty in the mathematical description of this physical motion is due to the origin of the field, which consists of the contribution of mobile ions and of a possible external field. The dielectric membrane also affects the field by image charges. So far only few electro-diffusion systems are well understood: although the voltaic cell was invented more than 200 years ago, designing optimal configurations is still a challenge. On the other extreme, ionic flux and gating of voltage-channels [2] is now well explained by the modern Poisson-Nernst-Planck theory of electro-diffusion [16], because at the nanometer scale, the cylindrical geometry approximation of protein channels reduces the computation of the electric field and of ionic diffusion to one dimension [9, 21, 15, 8, 7, 24, 25]. However, cellular domains at a micron scale involve two- and three-dimensional geometry, much more complicated than the cylindrical geometry of a channel pore, leading to a more complex electro-diffusion description [23, 14].

We recall that local curvature is a key geometrical element for controlling charge distribution in various media, such as in the air (e.g., the lightning rod [5]). The manifestation of this effect is observed in Lebesgue’s thorn, which is a an inverted cusp singularity of the boundary, for which the solution of Laplace’s equation blows-up inside the domain [5, p.304]. In electronics, the design of printed circuits is always pre-conditioned on corner effects [22]. However, these effects are not very well known inside an electrolytic bath. Recent analysis [4], [14] suggests that non-electro-neutrality in the geometry of an electrolyte confined by a dielectric membrane affects charge distribution.

We use the Poisson-Nernst-Planck (PNP) equations for charge concentration and electric potential as a model of electro-diffusion. The entire boundary is impermeable to particles (ions) and the electric field satisfies the compatibility condition of Poisson’s equation. Phenomenological descriptions of electro-diffusion, such as the cable equation or the reduced electrical engineering approximation by resistance, capacitance, and even electronic devices, are not sufficient to describe non-cylindrical geometry [14], because they assume a simple reduced one-dimensional or reduced geometry. We present here results about charge and field distributions in electro-diffusion in various geometrical microdomains, when the condition of electro-neutrality is not satisfied. We recall that under the non-electro-neutrality assumption, and with charge distributed in bounded domains confined by a dielectric membrane, Debye’s concept of charge screening decaying exponentially away from a charge [6], do not apply and long-range correlation leads to a gradient of charges in a ball with no inward current. A new capacitance law was derived for an electrolyte ball [4], where the difference of potential between the center $C$ and the surface $S$, that is, $V(C)-V(S)$, increases, first linearly and then logarithmically, when the total number of charges in the ball increases.

Our aim here is to understand the effect of boundary curvature on an electrical cell, such as neuron. In particular, we explore the effect of boundary curvature on the charge and field distribution at steady state. The curvature of membranes of dendrites and axons of neurons have many local maxima that can modulate the channel’s local electric potential [27]. In this article, we study the effects of local curvature on the distribution of charge in bounded domains with no electro-neutrality. The effect of non-electro-neutrality was recently studied in [4] and a long-range electrostatic length, much longer than the Debye length was found. This effect is due to the combined effects of non-electro-neutrality and boundary, which lead to charge accumulation near the boundary.

The cusp-shaped funnel geometry was studied in [10], however this paper presents several crucial mathematical differences with [10], in particular, we are solving a nonlinear equation, while it was linear in [10]. Furthermore, the boundary condition at the end of the cusp-shaped funnel: while it is the Dirichlet condition in [10], it is the Neumann condition here. This means that in [10] the absorption flux at the end of the funnel is computed, whereas here the stationary voltage and charge distribution are computed in the absence of flux. We develop here new boundary layer analysis, different than the classical matched asymptotics method [28, 29, 30]. The manuscript is organized as follow: first, we consider a bounded domain with an uncharged narrow cusp-shaped funnel on the boundary, which is a singular geometrical effect. Second, we further study the case of charge distribution in a charged narrow cusp.

2 The PNP equations

The Poisson-Nernst-Planck system of equations in a domain $\Omega$, whose dielectric boundary $\partial\Omega$ is represented as the compatibility condition for Poisson’s equation, and its impermeability to the passage of ions is represented as a no-flux boundary condition for the Nernst-Planck equation. We assume that the total charge in $\Omega$ consists of $N$ identical positive ions with initial particle density q(\mbox{\boldmathx}) in $\Omega$, their valence is $z$, and the total number of particles is fixed, equal to

[TABLE]

Thus the charge in $\Omega$ is

[TABLE]

where $e$ is the electronic charge. The charge density \rho(\mbox{\boldmathx},t) is the solution of the initial and boundary value problem for the Nernst-Planck equation

[TABLE]

Here \phi(\mbox{\boldmathx},t) is the electric potential in $\Omega$ and is the solution of the Neumann problem for the Poisson equation

[TABLE]

where \sigma(\mbox{\boldmathx},t) is the surface charge density on the boundary $\partial\Omega$. In the steady state,

[TABLE]

3 Steady solution in a ball with a cusp-shaped funnel

Local boundary curvature is a key geometrical feature that controls charge distribution in the domain. Specifically, we study the effect of a narrow funnel attached to a sphere. In various media, such as air (e.g., the lightning rod, [5]), the manifestation of this effect is observed in Lebesgue’s thorn, which is a an inverted cusp singularity of the boundary, for which the solution of Laplace’s equation blows-up inside the domain [5, p.304]. In the steady state (2) gives the particle density

[TABLE]

hence (5) gives Poisson equation

[TABLE]

and (6) gives the boundary condition

[TABLE]

for |\mbox{\boldmathx}|=R, which is the compatibility condition, obtained by integrating Poisson’s equation (5) over $\Omega$. Changing variables to

[TABLE]

Poisson’s equation (9) becomes

[TABLE]

and the boundary condition (10) becomes

[TABLE]

The translation \tilde{u}=u+{\displaystyle\ln\left(\lambda/\int_{\Omega}\exp\{v(\mbox{\boldmathx})\}\,d\mbox{\boldmathx}\right)}, converts (12) into

[TABLE]

We consider a dimensionless planar domain $\Omega$ with a cusp-shaped funnel formed by two bounding circles $A$ and $B$ of dimensionless radii 1 (see Fig.1(left)). The opening of the funnel is $\varepsilon\ll 1$. We construct an asymptotic solution in this limit to the nonlinear boundary value problem (BVP) (14) by first mapping the domain $\Omega$ conformally with the Möbius transformation of the two osculating circles $A$ and $B$ into concentric circles (see Fig.1(right)). To this end, we move the origin of the complex plane to the center of the osculating circle $B$ and set

[TABLE]

where

[TABLE]

The Möbius transformation (15) maps the circle $B$ (dashed blue) into itself and $\Omega$ is mapped onto the domain $\Omega_{w}=w(\Omega)$ in Figure 1(right).

The straits in Figure 1(left) are mapped onto the ring enclosed between the like-style arcs and the large disk is mapped onto the small red disk in Figure 1(right). The radius of the small disk and the elevation of its center above the real axis are $O(\sqrt{\varepsilon})$. The short black segment $AB$ of length $\varepsilon$ in Figure 1(left) is mapped onto the segment $AB$ of length $2\sqrt{\varepsilon}+O(\varepsilon)$ in Figure 1(right). This mapping (see [11]), transforms the PNP equations as well and thus leads to a new non-linear effect. Setting $u(z)=v(w)$ converts (12) to

[TABLE]

The boundary segment $AB$ at the end of the cusp-shaped funnel in Figure 1(left) is denoted $\partial\Omega_{w,a}$. To determine the boundary conditions, we use the change of coordinates $w=Re^{i\theta}=X+iY$. At the end of the funnel, where $R\simeq 1$, we get

[TABLE]

where

[TABLE]

For $\theta=\pi$ (for $z=-1$), we obtain $\partial\theta/\partial Y=-2/\sqrt{\varepsilon}$ and the boundary condition at $\partial\Omega_{w,a}$ is

[TABLE]

3.1 Reduced PNP equations in an uncharged cusp-shaped funnel

Approximating the banana-shaped domain $\Omega_{w}$ by a one-dimensional circular arc, we use a one-dimensional approximation of the solution in $\Omega_{w}$ [12, 13]. This approximation assumes that there are no non-neutralized charges on the surface of the cusp (Fig.3A). The boundary condition for the approximate one-dimensional solution of (17) is zero at angle $\theta_{Lim}=c\sqrt{\varepsilon}$, where $c$ is a constant (see details in [12, 13]) and represents the solution inside the disk in Figure 1(left), away from the cusp. Thus, (17) in the conformal image $\Omega_{w}$ becomes the boundary value problem

[TABLE]

The solution of (21) is shown in Figure 3B-C in the two domains, $\Omega$ (panel A) and its image $\Omega_{w}$ (panel B).

Our goal is now to estimate the difference of potentials between the north pole $N$ and the end of the funnel $C$,

[TABLE]

To construct an asymptotic approximation to the solution of (22) in the limits $\varepsilon\to 0$ and $\lambda\to\infty$, we first construct the outer-solution in the form of a series in powers of $\varepsilon$, which is an approximation valid away from the boundary. In the limit of small $\varepsilon$, the first term in the series vanishes, exponential terms drop out, and the second order term is

[TABLE]

where $M$ and $M^{\prime}$ are yet undetermined constants. The outer solution cannot satisfy all boundary conditions, so a boundary layer correction is needed at the reflecting boundary at $\theta=c\sqrt{\varepsilon}$. Thus, we set $\theta=\sqrt{\varepsilon}\xi$ and expand

[TABLE]

Writing the boundary layer solution as $y_{\mbox{\scriptsize bl}}(\theta)=Y(\xi)$, we obtain to leading order the boundary layer equation

[TABLE]

with $Y^{\prime}(c)=0$. The solution is decaying for large $\xi$ and develops a singularity at finite $\xi$. However, a Taylor expansion near $\xi=0$,

[TABLE]

gives in (25)

[TABLE]

In general, the coefficients satisfy $B_{k}=O(e^{-A})$, for $A\gg 1$. For small $\xi$, we obtain the approximate solution of (25) by considering the leading term in a regular expansion of the solution in powers of $\xi$. The equation for the leading term is

[TABLE]

and the solution is defined up to an additive constant. Setting $Y_{appr}(0)=0$, which does not affect the potential difference, we find that

[TABLE]

It follows that the boundary layer solution at $c\sqrt{\varepsilon}$ is

[TABLE]

The boundary layer near $\pi$ is needed, because $A\to\infty$ as $\varepsilon\to 0$ (see (45) below). An approximation of the solution can be obtained by freezing the power-law term in (22), for which the equation is for a generic parameter $b>0$,

[TABLE]

The solution is

[TABLE]

Putting the outer and boundary layer solutions together gives the uniform asymptotic approximation

[TABLE]

where the parameters $A$ and $b$ are yet undetermined constants. The condition at $c\sqrt{\varepsilon}=o(1)$ for $\varepsilon\ll 1$ is satisfied, because

[TABLE]

The condition at $\theta=\pi$ gives that

[TABLE]

The compatibility condition for (14),

[TABLE]

gives in $\Omega_{w}$ that

[TABLE]

Using the uniform approximation (32) in the compatibility condition (34), we obtain the second condition

[TABLE]

where we used the change of variable $\theta=\sqrt{\varepsilon}\xi$. Integrating by parts, we get for $\varepsilon\ll 1$

[TABLE]

where

[TABLE]

Thus, it remains to solve the asymptotic equation

[TABLE]

for $A$ and $b$ in the limit $\varepsilon\to 0$. We consider the limiting case where

[TABLE]

for which condition (3.1) can be simplified and gives to leading order

[TABLE]

that is, for $\lambda\sqrt{\varepsilon}\ll 1$ (40) gives

[TABLE]

With condition (38), we get

[TABLE]

To leading order in large $C$, we obtain

[TABLE]

The solution is expressed in terms of the Lambert-W function,

[TABLE]

and for small $\varepsilon$, using the asymptotics of the Lambert function,

[TABLE]

Finally,

[TABLE]

It follows that a uniform asymptotic approximation (32) in the limits $\lambda\to\infty$ $\varepsilon\to 0$ is given by

[TABLE]

The uniform approximation (46) is plotted for different values of $\varepsilon$ and $\lambda$ in Figure 3 against the numerical solution of (21), with the boundary conditions $v^{\prime}(c\sqrt{\varepsilon})=v^{\prime}(0)=0$. The numerical solutions are computed with the software COMSOL, based on an adaptive mesh refinement and a relative tolerance of $10^{-3}$, that we validated on known analytical results of steady state PNP equations in a disk [4]. We find that the asymptotic expansion is particularly good in the limit $\varepsilon\rightarrow 0$ and $\lambda\rightarrow\infty$ (Fig.3A-D). However, for $\lambda=O(1)$ the log-term approximation in (46) is non-monotonic in $\theta$.

Finally, to further validate the uniform asymptotic expansion, we compared the numerical solutions of the full equation (23) in the initial domain $\Omega$ with the reduced PNP equation (14) with zero Neumann boundary conditions, except at the end of the funnel for the mapped domain $\Omega_{w}$. The result is shown in Figure 3A-D, showing good agreement between the one-dimensional PNP approximation in $\Omega_{w}$ and the numerical solution of the full equation. We estimated numerically the difference between the asymptotic solution $y_{\mbox{\scriptsize unif}}$ (32) and the numerical estimation $V_{num}$ (21), averaged over the domain $\Omega_{w}$, for $10^{3}\leq\lambda\leq 5\cdot 10^{4}$ and $5\cdot 10^{-3}\leq\varepsilon\leq 10^{-1}$. The difference is almost constant in the range $[0.01,0.025]$.

To compare the voltage at the north and south poles (at the end of the funnel), we use the two-dimensional analytical solution in the entire ball and the numerical solution of (14) (Fig.3D). Interestingly, we find that the difference $u(N)-u(S)$ has a maximum with respect to $\lambda$, where $u(N)$ and $u(S)$ are the values of the potential at the north pole and at the end of the funnel, respectively.

3.2 The voltage drop between the end of the funnel and the center of the ball

We can now use (32) to compute the potential drop in (23). It is given by

[TABLE]

The potential difference $\tilde{\Delta}_{SC}u$ with respect to $\lambda$ is shown in Figure 3F (red line).

Next, we compare the potential drop (23) with the one between the center and the north pole. Numerical solution of the PNP equations shows that the voltage and charge distribution in a disk with a funnel do not differ from the ones in a disk in the upper sphere (Fig.2). This result is compared next to the difference between the north pole and the center evaluated from the exact analytical expression derived for a disk.

The expression for the voltage in the two-dimensional disk of radius $R$ is given by (see [4])

[TABLE]

where $\lambda_{D}$ is a parameter. We calibrate $\lambda_{D}$ so that the solutions of the PNP equations in a disk with a funnel have the same total charge as a disk. The Neumann boundary conditions for the disk and the funnel are, respectively,

[TABLE]

The calibration is

[TABLE]

We compare in Figure 3D the two-dimensional numerical solution of the PNP equation (14) in the domain $\Omega$ (blue line), with the analytical solution (3.2) in a disk with no cusp (dashed red). The numerical solution of the PNP equation (14) is plotted along the main axis $0y$ in the interval $[0,y_{0}]$ (where the point $y_{0}$ is defined by the condition $\nabla u(y_{0})=0$). In the range $[y_{0},y_{cusp}]$, where $y_{cusp}$ is the coordinate of the cusp, we compare the solution of (14) with the uniform solution $y_{unif}$ of (32) in the funnel (dashed green). We conclude that in the cusp, the two-dimensional approximation in a disk is in good agreement with the numerical solution of equation (14), confirming that the solution in the bulky head does not influence the one in the cusp (as already shown in Fig. 2). This result also confirms the validity of the analytical formula to predict the large $\lambda$ asymptotics.

For a disk of radius $R$, the potential drop is given by

[TABLE]

(see section 3.2). The potential drop $\tilde{\Delta}_{NC}u$ is shown in Figure 3E (blue line). The two differences of potential $\bar{\Delta}_{SC}u$ (49) and $\tilde{\Delta}_{NC}u$ (3.2) have the same logarithmic behavior $\ln{1/\lambda^{2}}$ for $\lambda\gg 1$ and $u(N)-u(S)=O(1)$. A numerical solution in two-dimensions shows that $u(N)-u(S)$ may converge to zero as $\lambda$ increases (Fig. 3F), thus having a local maximum for small values of $\lambda$. This maximum cannot be analyzed by the uniform expression (32), because it appears outside the domain of validity of (32). This result is in agreement with the two-dimensional numerical solution of (14) for the difference between $u(N)$ (potential at the north pole) and $u(S)$ (potential at the end of the funnel) (Fig. 3F).

The potential drop calculated above is non-dimensionalized by the radius of curvatures $R_{f}$ at the right and left of the funnel,

[TABLE]

where $\tilde{\varepsilon}$ is the length of the absorbing arc $AB$. The non-dimensionalized volume and boundary measure are, respectively,

[TABLE]

In dimensional units (3.2) gives the potential drop in the dimensional disk with a funnel as

[TABLE]

We conclude in the limit of $\lambda\gg 1$, $\tilde{\varepsilon}\to 0$ that the difference of potential between the end of cusp $S$ and the north pole $N$ in the domain is obtained by adding (49) and (50) and we get

[TABLE]

We recall that $R$ is the radius of the entire ball, while $R_{f}$ is the radius of curvature of the funnel.

4 The PNP equations in a charged domain with a cusp-shaped funnel

Due to the Neumann boundary condition (6) on the lateral part of the funnel,(17) in the transformed domain cannot be reduced to one dimension. Thus we derive a different one-dimensional approximation for the mapped PNP equations in the banana-shaped domain $\Omega_{w}$ by averaging over the radius $r$. Rewriting (17) in polar coordinates $w=re^{i\theta}$, we obtain

[TABLE]

In the section $\Omega_{w}\cap\{1-\sqrt{2\varepsilon}<r<1]\}$, the boundary conditions are

[TABLE]

Taylor’s expansion of $v$ in the section gives

[TABLE]

and because $|r-1|=O(\sqrt{\varepsilon})$, we obtain the approximation,

[TABLE]

Multiplying (52) by $r^{2}$ and integrating over the radius, we get

[TABLE]

The boundary conditions (53) give, to leading order in $\sqrt{\varepsilon}$, that

[TABLE]

that is, the BVP (52) in the section becomes the ODE (with respect to $\theta$),

[TABLE]

The graph of the solution of (57) in $\Omega_{\omega}$ and $\Omega$ is shown in Figure 6. Equation (57) is obtained by averaging over the radial direction and its solution seems to be a good approximation to (52) only for small $\lambda$. A different approach for large $\lambda$ is discussed in the next section.

A regular expansion for $\lambda\ll 1$,

[TABLE]

gives in (57) that $w_{0}=O(\varepsilon)$ and $w_{1}$ is the solution of the BVP

[TABLE]

Direct integration with respect to $\theta$ gives

[TABLE]

Equation (60) gives $A$ as

[TABLE]

The zero Neumann boundary condition cannot be satisfied and a boundary layer appears, leading to the local expansion

[TABLE]

It follows that for $\lambda\ll 1$, the solution increases with $\lambda$. It is shown below that it decreases for $\lambda\gg 1$, demonstrating that there is at least one maximum in the variable $\lambda$.

4.1 PNP asymptotics in a charged disk with a charged funnel

In the limit of $\lambda\gg 1$, $\varepsilon\to 0$, the asymptotic expansion of the potential found above for a charged disk with a funnel is no longer valid. Some insight can be gained by observing the field lines in the domain $\Omega_{w}$, described in Figure 2B. These lines are parallel to the radius vector, except in a small region near the funnel. Two sections can be distinguished,

[TABLE]

The two sections $A$ and $B$ are illustrated in Figure 5A. Note that the boundary of section $B$ contains a circular arc (marked magenta). Next, the approximate solutions $u_{A}(r,\theta)$ and $u_{B}(\theta)$ of (52) in the two sections are constructed

and used to construct a uniform approximation $u_{unif}$ in $\Omega_{w}$ (Fig. 5B).

4.2 Asymptotics of $u_{A}(r,\theta)$ in section $A$

The boundary conditions (53) for the potential equation (52) indicate that the radial derivative is $O(\lambda\sqrt{\varepsilon})\to\infty$ Thus the angular derivatives are negligible relative to the radial ones. It follows in a regular expansion of the solution that the $\theta$ derivatives can be neglected relative to the $r$ derivative and the equation is then solved along the rays $\theta=const=\theta_{0}$ for $r\in[1-\sqrt{\varepsilon},1]$. Thus, to leading order in $\lambda\sqrt{\varepsilon}$,

[TABLE]

For $\varepsilon\ll 1$, we get $|re^{i\theta_{0}}(1-\sqrt{\varepsilon})-1|^{4}=|e^{i\theta_{0}}-1|^{4}+O(\sqrt{\varepsilon})$. Setting

[TABLE]

and

[TABLE]

we get

[TABLE]

The general solution of (68) is given by [4]

[TABLE]

where the constants $C_{1}$ and $C_{2}$ are determined from the boundary conditions (68). Using

[TABLE]

we find the constant $C_{1}$ from (70) and from the boundary condition (68) at the point $r=1-\sqrt{\varepsilon}$, getting

[TABLE]

This gives in (70) at $r=1$ the transcendental equation for $C_{2}$,

[TABLE]

hence

[TABLE]

Now, it follows from (71) that

[TABLE]

Note that $\lim_{\lambda\rightarrow\infty}C_{2}\neq 0$, because otherwise we would get the asymptotic expansion

[TABLE]

which leads to

[TABLE]

and contradicts the condition (72) in the limit $\lambda\to\infty$.

Then (73) and (74) would imply that

[TABLE]

and (77) would give $C_{2}\gg 1$, so that the $\arctan$ term in (74) drops out, and we would be left with

[TABLE]

hence

[TABLE]

Expanding the left hand side of (72), using (73) and (78), we obtain that

[TABLE]

Together with (80), the solution of (72) is

[TABLE]

With the values of $C_{1}$ and $C_{2}$ computed in (71) and (81), the solution $v_{A,\theta_{0}}$ of (69) is given by

[TABLE]

Finally, using (67) and (82), we obtain for ($r,\theta)\in A$,

[TABLE]

The asymptotic solution $u_{A}$ is plotted in Figure 5B (blue dashed line). Comparison with numerical solutions for various values of $\lambda$ and $\varepsilon$ is shown in Figure 6 below.

4.3 The asymptotics of $u_{B}$ in section $B$

The asymptotic solution $u_{A}(r,\theta)$ in section $A$ cannot satisfy the boundary condition (53) at $\theta=\pi$. Indeed, (83) gives $\partial u_{A}(r,\theta)/\partial\theta|_{\theta=\pi}=0$, while the boundary condition (57) is $\partial v/\partial\theta|_{\theta=\pi}=-\lambda\sqrt{\varepsilon}/2|\partial\Omega|$, so a boundary layer correction is needed.

The boundary layer $u_{B}(\theta)$ is an asymptotic solution of (52) in section $B$, where the $\theta$ derivatives dominate the radial ones. The right-hand-side of (52) can be simplified for $\varepsilon\ll 1$.

For $r=1-\sqrt{\varepsilon}$ the approximation

[TABLE]

holds, which does not depend on $r$ and $\theta$. With this simplification in (52), we rewrite $u_{B}(\theta)$ as

[TABLE]

where $C_{0}$ is an additive constant and $\tilde{u}_{B}$ is a function of $\eta=\theta-(\pi-\sqrt{\varepsilon})$ and solves the BVP

[TABLE]

The solution of (86) (see [4]) is

[TABLE]

where $I_{\lambda}$ is the solution of the transcendental equation

[TABLE]

We obtain to leading order for $\lambda\gg 1$ that

[TABLE]

It follows from (89), (87), and (85) that for $\theta\in B$, the asymptotic solution is

[TABLE]

(see (85)). It is shown in Figure 5B (red dots).

4.4 A uniform approximation of $u(r,\theta)$ in $\Omega_{w}$

A uniform asymptotic approximation $u_{unif}(r,\theta)$ of the voltage $u(r,\theta)$ in the entire mapped domain $\Omega_{w}$ can be now constructed by matching the the leading term $u_{A}(r,\theta)$, given in (83) in section $A$, with that of $u_{B}(\theta)$, given in (90) in section $B$.

These approximations agree at $\theta=\pi-\sqrt{\varepsilon}$, so we obtain that

[TABLE]

Thus

[TABLE]

The numerical solution of (14) in $\Omega_{w}$ and the approximation $u_{unif}(r,\theta)$ of (92) are shown Fig. 5C.

4.5 Potential drop in $\Omega_{w}$

The potential drop $\tilde{\Delta}_{funnel}u$ between the center of mass $C$ and the tip of the funnel $S$, is

[TABLE]

Due to the axial symmetry of the domain $\Omega$, the center of mass $C$ is at $r=1-\sqrt{\varepsilon}$, hence (92) gives

[TABLE]

Recall that the constant $c$ depends on the domain geometry only, and is defined by the conformal mapping $w$ (see relation (15)). The potential drop $\tilde{\Delta}_{Cusp}u$ in the funnel can be decomposed as the sum of difference of potential between the two sections, $A$ and $B$. First, the approximations are

[TABLE]

and

[TABLE]

so that

[TABLE]

Using (83) in $A$, we get

[TABLE]

For $\varepsilon\ll 1$, we get from (4.5) that

[TABLE]

Hence, using (99) in (4.5), we get

[TABLE]

The approximate solution $u_{A}(S)$ at the tip of the funnel $S$ (south pole at $\theta_{0}=\pi$) is (100)

[TABLE]

At the center $C$, where $\theta_{0}=c\sqrt{\varepsilon}$, equation (4.5) gives for $\varepsilon\ll 1$ the $\theta_{0}$-dependent terms in (100) as

[TABLE]

and

[TABLE]

Using (102) and (103), the expression (4.5) reduces to

[TABLE]

For $\lambda\gg 1$, (104) becomes

[TABLE]

Finally, the approximate potential difference $\tilde{\Delta}u_{A}$ in (95), is the difference between (105) and (101),

[TABLE]

For $\lambda\gg 1$ (106) becomes to leading order

[TABLE]

which is independent of $\lambda$. (90) shows that the approximate potential in section $B$ is

[TABLE]

and

[TABLE]

Using (108) and (109) in (96), we obtain

[TABLE]

For $\lambda\gg 1$, (110) shows that $\tilde{\Delta}u_{B}$ is

[TABLE]

Finally, using (106), (110) and (97), we find that the potential drop is

[TABLE]

Again, using (107), (111) and (97) for $\lambda\gg 1$ limit, we get the approximate potential drop as

[TABLE]

Equation (110) shows that for $\lambda\gg 1$, the potential drop in the funnel domain occurs mostly in the region $B$. The expression (4.5) is plotted in Figure 5A-D (red) and compared to $\ln\lambda^{2}+const$ (green) and to a two-dimensional numerical solution. The good agreement confirms the validity of the asymptotic expansion and thus confirming the new asymptotic formulas derived here. We conclude with the general formula for a dimensional cusp-shaped funnel where $|\partial\Omega|=\frac{|\partial\tilde{\Omega}|}{R_{c}}$ and $R_{c}$ is the radius of curvature at the cusp

[TABLE]

4.6 Expansion of the potential drop between $N$ and $S$

To expand the potential difference $u(N)-u(S)$ between the funnel tip $S$ and the north pole $N$ of $\Omega$, we first use the results (114) computed above, to expand the difference $u(C)-u(S)$, and then subtract (114) and (49). The the terms $2\ln(\lambda)$ drop out and we have

[TABLE]

where $R$ is the distance between the north pole $N$ and the center of mass $C$ and $R_{c}$ is the radius of curvature at the cusp. We obtain to leading order

[TABLE]

which is a constant that depend only on the center of mass $C$.

5 Discussion and conclusion

We have derived here new electrostatic laws in non-neutral confined electrolytes from nonlinear electro-diffusion theory (PNP equations). The effect of local geometrical structure, such as the local curvature of the boundary emerges from the asymptotic solution of the model. The PNP equations describe the charge concentration and electric potential. The new electrical laws are derived in the context of non-electro-neutrality, where we use a single ionic species. The approximation of the steady-state solution in a ball with an attached cusp-shaped funnel on its surface is new and the construction of the asymptotic expansion uses a new boundary layer analysis.

Using asymptotic and numerical solution of the PNP equation, we found here that for a sufficiently high number of charges, the charge concentration peaks at the end of the funnel in a charged funnel boundary domain, but this is not the case for an uncharged funnel domain (Fig.7A-C). This effect is clearly the result of the cusp-shaped geometry. The present analysis reveals that the curvature affects the membrane potential. We also find that the voltage increases logarithmically in the total number of excess charges $N$, which is valid for uncharged (3.2) and charged (52) cusp-shaped funnel on the boundary. We studied here the voltage changes and electro-diffusion under an excess of positive ions. The voltage difference in the limit $\lambda\rightarrow\infty$ is probably attenuated in a mixed ionic solution, but the electro-neutrality remains broken. Cytoplasmic ions are characterized by the following concentrations Na+ = 148ml, K+ = 10ml and Cl- = 4ml. There is a clear unbalance toward positive charges, however there are probably molecules of various sizes with negative charges to re-balance the charges. However, the motility of these proteins should be driven by a diffusion coefficient smaller than the one of the ions. This difference of mobility is certainly a key feature in maintaining non-electro-neutrality and then tuning the value of $\lambda$. However, in a system containing an excess of positive and a small amount of negative charges, we show in Appendix that the limit of the PNP equations in the bulk, when the number of negative charges tends to zero, can be obtained by a regular expansion of the solution. This result shows that the small amount of negative charges does not perturb much the distribution of positive ones. red Finally, note that we did not consider here nanometer structures, such as ionic channels, where a negative ionic charge can affect the motion of the other ions in the channel pore.

We conclude that local geometrical properties, such as curvature, can modulate the local voltage in biological cellular electrolytes when electro-neutrality is violated. This result generalizes the case of a ball, where the distribution of charges accumulates on the surface as the total charge increases [4]. Following a non uniform boundary curvature, we expect that charges will be non-uniformly distributed, leading to a difference of potential across the membrane with charges on its surface. Since, this difference of potential plays a key role in information processing at synapses, we conclude that the spine geometry, in particular its curvature may impact the coding or decoding of voltage through current [27].This effect may as well influence the propagation and genesis of local depolarization [20, 19, 14]. More realistic funnels, with two different curvature radii can be incorporated to the formalism presented by modifying the parameter $\alpha$ (15) as shown in [12]. The formalism presented in this paper can be applied beyond physiology, in particular in the design of nanopipettes with an optimal shape [18, 12] by modulating $\alpha$ (15) or with a patterned surface [26] by changing the surface charge density via $\lambda$ in region $A$ (65).

6 Appendix

red

Regular expansion of the PNP solution when there are an excess of positive and a small number of negative charges

We show that for the concentrations of cations and anions found in literature [31], the leading order solution of the electrical potential in the bulk can be obtained by considering positive charges only. We assume that the charge of an electrolyte confined in $\tilde{\Omega}$ consists of identical $N_{p}$ positive and $N_{m}$ negative ions with density q_{p}(\mbox{\boldmathx}) and q_{m}(\mbox{\boldmathx}) such as

[TABLE]

where $p$ and $m$ are positive and negative species respectively. The total charges in $\tilde{\Omega}$ is the sum

[TABLE]

The associated charge densities \rho_{p}(\mbox{\boldmathx},t) and \rho_{m}(\mbox{\boldmathx},t) satisfy the boundary value problem for the Nernst-Planck equation

[TABLE]

where $z_{i}$ is the valence and $D_{i}$ is the diffusion coefficient for the ion specie $i$. The electric potential \phi(\tilde{}\mbox{\boldmathx},t) in $\tilde{\Omega}$ is solution of the Neumann problem for the Poisson equation

[TABLE]

where \tilde{\sigma}(\tilde{}\mbox{\boldmathx},t) is the surface charge density on the boundary $\partial\tilde{\Omega}$. At steady-state, (119) gives

[TABLE]

where $\rho_{i,0}$ is obtained from no-flux boundary condition (120), thus

[TABLE]

Using the non-dimensionalized potential \displaystyle\tilde{u}(\tilde{}\mbox{\boldmathx})=\frac{e\,\phi(\tilde{}\mbox{\boldmathx})}{k_{B}T}, equation (123) becomes

[TABLE]

Using (122) and (125), we obtain

[TABLE]

where $l_{B}$ is the Bjerrum length. Using \displaystyle\mbox{\boldmathx}=\frac{\tilde{}\mbox{\boldmathx}}{R_{c}} and $\tilde{u}(\tilde{x})=u(x)$ where $R_{c}$ is the cusp curvature radius, (126) becomes

[TABLE]

The small parameter is $\zeta=\displaystyle\frac{N_{m}}{N_{p}}\ll 1$ because in the bulk, the concentration of negative charges such as chloride (about $4$ mM) is much smaller than positive ones (potassium and sodium account together roughly for $167$ mM [31]). A regular expansion of u(\mbox{\boldmathx}) is

[TABLE]

Using (128) in (127), in small $\zeta$ limit, we have

[TABLE]

and in $\Omega$

[TABLE]

which admit a regular solution. This result shows that the limit of the PNP equation when $\zeta$ tends to zero (small charge limit) gives v_{0}(\mbox{\boldmathx}), and thus we conclude that a small amount of negative charges does not perturb the distribution of the total excess of positive charge in the bulk.

The numerical procedure

Numerical solutions were constructed by the Comsol Multiphysics 5.0 (BVP problems), Maple 2015 (Shooting problems) and Matlab R2015 (Conformal mapping). The boundary value problems in 1D, 2D, and 3D were solved by the finite elements method in the Comsol ’Mathematics’ package. We used an adaptive mesh refinement to ensure numerical convergence for large value of the parameter $\lambda$.

We solved the PDEs by the shooting procedure for boundary value problems using Runge-Kutta 4 method, as well as solvers from Maple packages.