Poisson-Nernst-Planck equations in a ball
aa r X i v : . [ phy s i c s . c l a ss - ph ] M a y Poisson-Nernst-Planck equations in a ball
Z. Schuss ∗ , J. Cartailler and D. Holcman † August 20, 2018
Abstract
The Poisson Nernst-Planck equations for charge concentration and electric po-tential in a ball is a model of electro-diffusion of ions in the head of a neuronaldendritic spine. We study the relaxation and the steady state when an initialcharge of ions is injected into the ball. The steady state equation is similar to theLiouville-Gelfand-Brat´u-type equation with the difference that the boundary con-dition is Neumann, not Dirichlet and there a minus sign in the exponent of theexponential term. The entire boundary is impermeable to the ions and the elec-tric field satisfies the compatibility condition of Poisson’s equation. We constructa steady radial solution and find that the potential is maximal in the center anddecreases toward the boundary. We study the limit of large charge in dimension 1,2and 3. For the case of a small absorbing window in the sphere, we find the escaperate of an ion from the steady density.
The non-linear system of Poisson-Nernst-Planck (PNP) equations has been widely usedto study properties of the electric field in local nanodomains such as ionic channels[21, 19, 10, 5, 6]. It was also used to simulate the equilibration of ions between largereservoirs through narrow necks [21, 11]. It is also possible to study the effect of interact-ing ions in ionic channels cite [27, 9].We use here PNP to study the distribution of charges at the micrometer scale level. In-deed, the stationary PNP equations with Neumann and no-flux conditions on the bound-ary of a finite domain Ω, respectively, describe the electrical potential and density ofcharge in Ω. They can be reduced to a Liouville-Gelfand-Brat´u-type equation for theelectric potential with however two major differences: first, the boundary condition on ∂ Ω is Neumann but not Dirichlet and second, there is a minus sign in the exponent,which is normalized over the domain Ω. We study here the solution of this equation in ∗ Department of Mathematics, Tel-Aviv University, Tel-Aviv 69978, Israel, [email protected] † Ecole Normale Sup´erieure, 46 rue d’Ulm 75005 Paris, France, [email protected]. ≤ λ . We construct asymptotic approximations of the solutions for small and large λ . Theone-dimensional case is solved explicitly and it is characterize by a log-singularity at theboundary that develop in the large λ limit. The explicit solution in two-dimension hasalso a singularity on the boundary. We also obtain a similar asymptotic behavior inthree-dimension, although the solution cannot be computed explicitly and we provide anasymptotic and numerical argument for large λ showing again a log-singularity at theboundary. We study the voltage change from the center and how it develops a boundarylayer for large λ . The voltage drop from the center to the boundary converges to a finitevalue as λ increases to infinity. We also apply the analysis of PNP to study idealized den-dritic spine structure, which consists of a spherical dielectric membrane filled with ionicsolution, connected to the dendrite by a cylindrical narrow neck. The paper is organizedas follow: in the first part we study asymptoticall and numericall PNP. In the second part,we estimate the current generate in an idealized spine (head connected by a cylinder).We derive the current generate by a spine. We show here how the head geometry controlsthe voltage, while the narrow neck radius control the current. We consider the Poisson-Nernst-Planck system in a ball Ω of radius R , whose dielectricboundary ∂ Ω is represented as the compatibility condition for Poisson’s equation and itsimpermeability to the passage of ions is represented as a no-flux boundary condition forthe Nernst-Planck equation. We assume that there are N positive ions of valence z in Ωand that there is an initial particle density q ( x ) in Ω such that Z Ω q ( x ) d x = N. (1)The charge in Ω is Q = zeN, where e is the electronic charge. The charge density ρ ( x , t ) is the solution of the Nernst-Planck equation D h ∆ ρ ( x , t ) + zekT ∇ ( ρ ( x , t ) ∇ φ ( x , t )) i = ∂ρ ( x , t ) ∂t for x ∈ Ω (2) D (cid:20) ∂ρ ( x , t ) ∂n + zekT ρ ( x , t ) ∂φ ( x , t ) ∂n (cid:21) = 0 for x ∈ ∂ Ω (3) ρ ( x ,
0) = q ( x ) for x ∈ Ω , (4)where the electric potential in Ω is φ ( x , t ) is the solution of the Poisson equation∆ φ ( x , t ) = − zeρ ( x , t ) εε for x ∈ Ω (5)2nd the boundary condition ∂φ ( x , t ) ∂n = − σ ( x , t ) for x ∈ ∂ Ω , (6)where σ ( x , t ) is the surface charge density on the boundary ∂ Ω. In the steady state andin spherical symmetry σ ( x , t ) = − Q πR . (7) In the steady state ∂ρ/∂t = 0 so (2) gives the density ρ ( x ) = N exp (cid:26) − zeφ ( x ) kT (cid:27)Z Ω exp (cid:26) − zeφ ( x ) kT (cid:27) d x , (8)hence (5) gives ∆ φ ( x ) = − zeN exp (cid:26) − zeφ ( x ) kT (cid:27) εε Z Ω exp (cid:26) − zeφ ( x ) kT (cid:27) d x . (9)In spherical symmetry in R d (9) can be written in spherical coordinates as φ ′′ ( r ) + d − r φ ′ ( r ) = − zeN exp (cid:26) − zeφ ( r ) kT (cid:27) S d εε Z R exp (cid:26) − zeφ ( rkT (cid:27) r d − dr < , (10)where S d is the surface area of the unit sphere in R d . The boundary conditions are ∂φ (0) ∂r = 0 , ∂φ ( R ) ∂r = − QS d R d − . (11)The inequality in (10) means that φ ( r ) has a maximum at the origin and decreases towardthe boundary (see Fig. 3A). We can normalize the radius by setting r = Rx for 0 < x < u ( x ) = zeφ ( r ) kT , λ = ( ze ) Nεε kT , (12)3o write (10) as u ′′ ( x ) + d − x u ′ ( x ) = − λ exp {− u ( x ) } S d R d − Z exp {− u ( x )) } x d − dx (13) u (0) = 0 , u ′ (0) = 0 . Incorporating the denominator of the RHS of (13) into the parameter λ by setting λ = µS d R d − Z exp {− u ( x ) } x d − dx, (14)we can write the initial value problem (13) as u ′′ ( x ) + d − x u ′ ( x ) = − µ exp {− u ( x ) } (15) u (0) = u ′ (0) = 0 . First, we show that solutions exist in dimensions 1 ≤ d ≤ µ in the range0 ≤ µ < µ ∗ for some positive µ ∗ . Solution in dimension one
We solved directly equation 15 in dimension 1 (see appendix 6.1) and we obtain (see eq.52)that u Dλ ( x ) = ln cos r λ I λ x ! , (16)where I λ is solution of the implicit equation I λ = 2 λ tan r λ I λ . (17)The graph of u Dλ ( x ) is shown Fig. 2A, while the one for λI λ versus λ is shown in Fig. 2B.We have 0 < µ ( λ ) = λI λ ≤ π and lim λ →∞ µ ( λ ) = π . The solution exists u Dλ for all λ > x = 1 when λ → ∞ . Solution in dimension two
In dimension 2, we obtain the solution in (appendix 6.2) u Dλ ( x ) = log (1 − λ I λ x ) . (18)4here I λ = π + 18 λµ ( λ ) = λI λ lim λ →∞ µ ( λ ) = 8 . The graph of u λ ( x ) is shown on Fig. 2C, while the one for λI λ is on Fig. 2D. u λ ( x ) =log(1 − λλ +8 π x ) , develop a log-singularity as λ → ∞ . Analysis in dimension three
The solution of the initial value problem (13) in dimension d = 3 can be directly computed.We show now that the solution exits for all λ , while there is a critical value µ ∗ , abovewhich, there is no regular solution. Contrary to dimensions one and two, the value of µ ∗ can only be estimated numerically. We first show using a phase-space analysis that thesolution of equation 15 is unique when it exists. However it is not possible to use thephase-space to study the singularity of the equation. To study the asymptotic explosionof the equation, we use an asymptotic argument. Finally, we will study the solutionnumerically.Next, we show that the problem (13) has a unique regular solution for all λ ≥
0, whenthe solution is finite. The proof of uniqueness of the solution follows the phase-spaceanalysis of (15). Indeed, using the change of variables s = − log r, u ( r ) = U ( s ) , v ( s ) = dU ( s ) ds , w = µe − s e − U ( s ) w ′ ( s ) = − w ( s ) − U ′ ( s ) w ( s ) = w ( s )[ − − v ( s )] , (19)which gives v ′ ( s ) = v ( s ) − w ( s ) , w ′ ( s ) = − w ( s )[2 + v ( s )] , (20)and can be written as dwdv = − w (2 + v ) v − w . (21)The phase space of (20) contains exactly two critical points. The origin is a saddle pointand its stable manifold has the tangent T of equation w = 3 v . The point P a = ( − , − u (0) = u ′ (0) = 0 for the solution of (15)impose lim s →∞ U ( s ) = u (0) = 0 and lim s →∞ U ′ ( s ) = − lim r → ru ′ ( r ) = v (0) = 0, hencethe constraints lim s →∞ v ( s ) = 0 , lim s →∞ w ( s ) = lim s →∞ µe − s e − U ( s ) = 0 . (22)5hus the trajectory of the solution of (15) in the first quadrant, which satisfies the con-straints (22), has to be on the separatrix that converges to the saddle point. Choosingany value U (0) gives µe − U (0) the value of v (0) = U ′ (0) has to be chosen on the separatrix.Therefore starting in the first quadrant, a trajectory of (20) converges to the saddle pointif and only if it starts on the separatrix with the tangent T . The stable branch at thesaddle point tends to infinity as s decreases toward 0. Indeed, the local expansion of (21)near the saddle point is w ( v ) = 3 v + 35 v − v + . . . , (23)which gives the phase portrait (Fig. 1). Finally, along the separatrix w ′ ( v ) >
0, exceptat the origin, showing that for an initial v (0), there is a unique solution. However, itis not possible from the phase-space to study singular solution. Indeed, as we shall see,when there is a singularity, because the it occur precisely at the initial value and thus theCauchy problem cannot even start. We conclude the problem (13) has a finite solution −4 −2 0 2−10−50510 v w Q P
Figure 1:
Phase-space solution of (20). The separatrix is shown in red, while the othertrajectories are in blue. v (0) , w (0))(when it is finite) on the separatrix in the first quadrant, there is a unique solution to(20) that satisfies (22).A numerical solution of (13) gives the graph Fig. 3E, that is the solution u ( x ) of (3)for µ ≤ µ ∗ = 11 . µ ∗ = 14) blows up before reaching x = 1,while the dash (small point) graph is finite throughout the interval. To estimate an upperbound for µ ∗ , we note that whenever the solution exists for some µ near µ ∗ , its asymptoticbehavior for x close to 1 shows that u ′′ (1) ≫ u ′ (1) (see the blue graph in Fig. 3). Indeed,to show that under the assumption u ′′ (1) ≫ u ′ (1) the latter inequality is self-consistent,we note that near x = 1 the solution of (13) can be approximated by the solution of thesimpler problem ˜ u ′′ ( x ) = − µ exp {− ˜ u ( x ) } , (24)given by ˜ u ( x ) ∼ log cos (cid:18)r µ x (cid:19) . (25)Thus ˜ u ( x ) is finite in the interval as long as µ < π . µ ∗ (26)and ˜ u ′ ( x )˜ u ′′ ( x ) ≤ |√ µ − √ µ ∗ |√ µ ∗ ≪ . (27)We conclude at this stage that for fixed values of µ below and above, where above thelatter they blow-up inside the interval 0 < x < µ varies with λ according to (14), the solutions exist for all values of λ (frames B,D,F ofFig. 2). Figure 2A-C-E shows that potential drop between the center and the surface ofthe sphere as a function of λ for 1 ≤ d ≤ λ ≪
1, the ex-pansion u ( x ) = − λ x π + O ( λ ) (see eq. 74). In contrast, for λ ≫
1, we mention abovethat the approximation u ( x ) ≈ − x ), which was relevant near x = 1 can be usedin the entire domain [0 , The potential differences
The difference u (0) − u (1) as we shall see in the next section has a physical meaning, asit represents the difference of potential between the center and the periphery of a sphere.7 x u ( x ) µ = 2 µ = 4.5 µ = 6 A Dimension 1 µ lim = π /2 0 200 400024 λ λ / I λ B Dimension 1 π /20 0.5 1−10−50 x u ( x ) µ = 5 µ = 8 µ = 12 C Dimension 2 µ lim = 8 0 500 1000246 λ λ / I λ D Dimension 280 0.5 1−10−50 x u ( x ) µ = 8 µ = 10.5 µ = 14 E µ lim ≈ λ λ / I λ F Dimension 310.7
Figure 2:
Numerical solutions u ( x ) of the initial value problem (13). (A),(C), and (E)correspond to different profile values of λ in dimensions 1,2, and 3, respectively. The dottedcurves are solutions that blow-up for x <
1. (B),(D), and (F) are plots of the ratio λI λ indimensions 1,2 and 3, respectively. x u ( x ) λ =100 λ =1000 λ =10000 A −3 x AsymptoticsNumerics B λ << 1 0 0.5 1−15−10−50 x )Numerics C λ >> 1 Figure 3:
Asymptotics behavior of the solution u ( x ) . ( A) change in the profile u ( x ) for 3 values of the parameter λ = 10 , , . ( B,C)
We present two regimes: for λ = 0 . ≪
1, we have u ( x ) = − λ x π + O ( λ ) (see eq. 74) and λ ≫ u ( x ) ≈ − x ). The analytical approximations (red) are compared with the numerical solutions (seeappendix).We have in dimension 1, | u λ (1) − u λ (0) | = ln cos r λ I λ ! , (28)where λ I λ → π as λ → ∞ . in dimension 2, | u λ (1) − u λ (0) | = 2 log( 8 πλ + 8 π ) , (29)in dimension 3, for λ ≫ | u λ (1) − u λ (0) | = 2 log ln(1 − f ( λ )) , (30)where the function f is increasing and f ( λ ) → λ → ∞ . The different curves fordimension 1,2 and 3 are shown in fig. 4. In all cases, the large λ asymptotic is dominatedby the log-behavior. The distribution of voltage and charge in a dielectric ball can be estimated from the resultsof the previous results by using the dimensional relation 12 in a ball of radius R = 1 µm .We plotted in Fig. 5A the voltage drop for N = 10 , charges. Already for 1000charges, there is a difference between the center and the surface of a ball of few milli-Volt.This effect could be tested for in the head of dendritic spines. Moreover, the density ischarge is concentrated at the periphery (Fig. 5B), leading also to a large field close to the9
500 10000510 λ u ( ) − u ( ) Dim 1 Dim 2 Dim 3
Figure 4: Asymptotics of u λ (1) − u λ (0) for dimensions 1,2 and 3.boundary (Fig. 5C). Consequently most od the charge are accumulated at the boundary,as revealed by the plot of the cumulative density of charges Q ( r ) = N Z r exp (cid:26) − zeφ ( r ) kT (cid:27) r dr Z R exp (cid:26) − zeφ ( r ) kT (cid:27) r dr . (31)We conclude that when the total number of charges is fixed sufficiently high, the chargeaccumulate at the surface. The field is only significant close to the surface and thuscan trap a Brownian charge in such region, while outside a small boundary layer of theboundary, the field is almost zero and charge particle experience no drift. This effect isdiscussed in the section. Although we found previously that for a fixed radius, the difference of potential V (0) − V (1) is bounded as a function of the total number of charge, we shall now show thatthe maximal number of charges increases linearly with the radius of the ball. Indeed,introducing the dimensionless radial variable ζ = r/R and u λ ( r ) = U λ/R ( ζ ), equation (10)becomes U ′′ λ/R ( ζ ) + 2 ζ U ′ λ/R ( ζ ) = − λ exp (cid:8) − U λ/R ( ζ ) (cid:9) πR Z exp (cid:8) − U λ/R ( ζ ) (cid:9) ζ dζ , (32)10 Radius r ( µ m) V ( m V ) Charges10 Charges10 Charges A Radius r ( µ m) ρ ( µ m − ) Charges10 Charges10 Charges B Radius r ( µ m) E ( m V µ m − ) Charges10 Charges10 Charges C Radius r ( µ m) C u m u l . nb . c h a r g e s Charges10 Charges10 Charges D Figure 5:
Distribution of (A) the potential, (B) charge and the field (C) andcumulative density of charges (D) inside a dielectric ball. U λ/R (0) = U ′ λ/R (0) = 0. Now, we solve the initial value problem V ′′ µ ( ζ ) + 2 ζ V ′ µ ( ζ ) = − µ exp {− V µ ( ζ ) } , V µ (0) = V ′ µ (0) = 0 (33) W ′ µ ( ζ ) = ζ exp {− V µ ( ζ ) } , W µ (0) = 0and note that u λ ( r ) = V µ (cid:16) rR (cid:17) , λ = 4 πµRW (1) . (34)Thus the number of charges Q in a ball or radius R create the same distribution as acharge Q/R in a ball of radius one, which can be summarized as Q ( R ) = RQ (1) . (35) We now discuss various consequences of distributing charges close to the boundary, in thelarge charge regime. The first consequence is on the MFPT ¯ τ ( x ) from x ∈ Ω, which isthe solution of the Pontryagin-Andronov-Vitt (PAV) boundary value problem [20] D h ∆¯ τ ( x ) − zekT ∇ ¯ τ ( x ) · ∇ φ ( x ) i = − x ∈ Ω (36) ∂ ¯ τ ( x ) ∂n + zekT ¯ τ ( x ) ∂φ ( x ) ∂n =0 for x ∈ Ω r (37)¯ τ ( x ) =0 for x ∈ Ω a . (38)We consider the case of a large field −∇ φ ( x ) ≫ | x | = 1. The profileof φ ( x ) was studied in section 2.1 (see Figures 5). To study the solution of the PAVproblem (36)-(38), we map the neighborhood of ∂ Ω a smoothly into the upper half planewith coordinates X = ( x, y, z ), where z = 0 is the image of the boundary, ˜ τ ( X ) = ¯ τ ( x ),and outside a boundary layer near ∂ Ω a V = ∂φ ( x ) ∂n (cid:12)(cid:12)(cid:12) | x | =1 = const, Φ( x, y ) = φ ( x ) (cid:12)(cid:12)(cid:12) | x | =1 = const, so that ∇ x,y Φ( x, y ) = 0 . The PAV system (36)-(38) is converted to˜ uτ zz ( X ) − zekT V ˜ τ z ( X ) + ∆ x,y ˜ τ ( X ) = − D , (39)12 regular expansion of ˜ u ( X ) for large V gives that to leading order ˜ τ ( X ) is a functionof ( x, y ) and setting T ( x, y ) = ˜ τ ( x, y, x,y T ( x, y ) = − D . (40)Thus the MFPT from x ∈ Ω to ∂ Ω a is the sum of the MFPT from x to ∂ Ω and theMFPT form ∂ Ω to ∂ Ω a on the surface ∂ Ω. The MFPT to ∂ Ω is negligible relative to thatto ∂ Ω a . This approximation means that to reach ∂ Ω a in a highly charged ball a charge isfirst transported by the field to the reflecting part ∂ Ω r of the sphere with overwhelmingprobability and then it finds ∂ Ω a by surface diffusion. A second consequence of the charge distributions is the control of spine current, inde-pendently of the voltage. Indeed, the solution T ( x, y ) of (40) is the MFPT of Brownianmotion on a sphere of radius R to an absorbing circle centered at the north-south axisnear the south pole, with small radius a = R sin δ . It is given by [25] T ( x, y ) = 2 R D log sin θ sin δ , (41)where D is the diffusion coefficient, θ is the angle between x and the north pole. Thus¯ τ ( x ) = T ( x, y ) . (42)The MFPT, averaged over the sphere with respect to a uniform distribution of x is givenby ¯ τ = 2 R (cid:18) log 1 δ + O (1) (cid:19) for δ ≪ . (43)The MFPT for N independent charges is¯ τ N = 2 R N (cid:18) log 1 δ + O (1) (cid:19) for δ ≪ . (44)It follows that the current through the small window is given by J = ze ¯ τ N = QD R (cid:18) log Ra + O (1) (cid:19) for a ≪ R. (45)We conclude that once a current enters into a dielectrics ball such as a spine head, theexcess of charges Q is first pushed toward the boundary and before moving by Brownian13otion to the spine neck. This result shows that the current in a spine head is governedby the spine geometry and a key parameter is the radius a of the neck. When there is aconservation of charge principle (no leak), the current through the dendritic shaft is thesame as the one exiting the spine. In that conditions, the spine neck length do not affector modulate the current. Determining the voltage drop between the membrane of the spine head and the dendritewhen a current is flowing from the head to the dendrite remains challenging because thecable theory cannot be applied in a system that cannot be approximated by a cable.The general scheme for modeling the electro-diffusion in the spine is the PNP model inthe head and a one-dimensional conduction of ions in the neck. The neck is considereda classical ionic conductor. Thus the steady-state PNP equations have to be solved inthe sphere with boundary conditions implied by the compatibility condition and the fluxthrough the neck is determined by the mean first passage time (MFPT) of ions from thehead to the neck, as discussed above. In the case of high charge Q the potential turns outto be practically flat throughout the ball with a sharp boundary layer with negative slopeat the boundary. Thus charge diffuses and is pushed strongly toward the membrane soionic motion is practically confined to motion on the surface. Due to spherical symmetry,the potential is constant on the boundary so ionic motion is free Brownian motion on asphere. At high charge ions interact through the ambient potential that is determined fromPoisson’s equation in the ball. Therefore they can be assumed independent free Brownianparticles. The MFPT ¯ τ of an ion to the small opening of the neck is determined from thetwo-dimensional NET theory(see previous section). Because the flux carried by a singleion is q/ ¯ τ , where q is the ionic charge, the number of ions in the spine head is N = Q/q and the MFPT ¯ τ N of any of the N ions is given by¯ τ N = ¯ τN . Thus the current through the neck is I = Q ¯ τ (46)and due to charge conservation, it is independent of the length of the neck. If we considerthe neck to be a parallel-plate capacitor carrying a steady current I , then the voltagedrop across the neck is simply V = RI , where R is the resistance of the neck, given by R = k B T L q nD , (47)14here k B is Boltzmann’s constant, T is absolute temperature, L is the length of theconductor, n is the number of ions in the neck, q is the charge of an ion, and D is thediffusion coefficient of the solution in the neck [20]. This model is valid as along as thevoltage is maintained in the spine head. We have studied here the solution of the PNP equations in a ball. We estimated the voltagein dendritic spines when the voltage in the spine head is maintained. A certain fractionof spines receive synaptic connections, essential for neuronal communication. Althoughtheir functions are still unclear, there are involved in regulating synaptic transmissionand plasticity [29, 18, 26, 8]. Interestingly, most of the excitatory connections occurs noton the dendrite but rather on spines and the reason is still not clear. The spine shape isquite intriguing, made of a head connected to the dendritic shaft by a cylinder. We foundhere that this geometry play a key role: the spine head geometry determines the drop ofpotential, while the current is defined by the diffusion on the surface and the mean time tofind the entrance of the neck in a two dimensional Brownian motion (see Narrow escapetime [22][13]). In the neck, under a voltage clamp condition, when a constant voltagedifference between the head and the neck is imposed, the voltage-current relation followa resistance law. Thus the spine geometry defines both the capacitance and resistance ingeometrical terms, a vision that complement previous classical studies [26, 17, 23].Finally, computing in the transient regime, the change in voltage drop between thespine head and the dendritic shaft, requires computing the time dependent PNP equa-tions. Another open question is to study the influence of the spine head geometry on thedistribution of charge. Computing the distribution of charges and the associated field innon-convex geometry is certainly the most challenging.
In this appendix, we first solve analytically the Liouville equation 13 in dimensions oneand two and in the second part, we describe the numerical methods to compute thesolution in dimension 3.
Liouville equation in the interval [0 1] is − u ′′ ( r ) = λ e − u R e − u dr (48)15ith initial conditions u (0) = 0 and u ′ (0) = 0 . (49)This is the classical Cauchy problem. After a direct integration we get with the initialconditions u ′ ( x ) = 2 λI λ ( e − u ( x ) − , (50)where I λ = Z e − u λ ( x ) dx. (51)A second integration gives u λ ( x ) = ln cos r λ I λ x. (52)Now we self-consistently calculate I λ = Z e − u λ ( x ) dx = Z dx cos q λ I λ = 1 q λ I λ tan r λ I λ . (53)Thus I λ > I λ = 2 λ tan r λ I λ . (54)The graph of λI λ versus λ is shown in Figure 2. We have lim λ →∞ λI λ = π , and specifically, y λ = q λ I λ = π − π λ + O ( λ ). The solution (52) is shown in Fig. 2 and is regular in theentire interval 0 < x < λ . The drop between the extreme points of theinterval is u λ (1) − u λ (0) = ln cos r λ I λ (55)and becomes infinite as the total charge increases indefinitely. The dimension 2 case can be transformed into the one dimensional case [14] using thechange of variable r = e − t ˜ u ( t ) = u ( r ) − t. − ˜ u tt = λI λ e − ˜ u ( t )+2 t (56)where I λ = 2 π R e − u ( r ) rdr and w ( t ) = ˜ u ( t ) + 2 t satisfies − w tt = λ e − w ( t ) I λ (57)The initial condition are now transform to asymptotic conditions at infinity:lim t →∞ ( w ( t ) − t ) = 0 (58)lim t →∞ ( ˙ w ( t ) − e t = 0 (59)A first integration gives ˙ w λ e − w ( t ) I λ + 2 . (60)The solution is w ( t ) = − log( 8( λe C +2 t − ) − C − t, (61)where C is a constant. Finally, we obtain that u λ ( r ) = log (1 − λ I λ r ) . (62)To close the equation, we shall now compute the integral I λ = Z e − u λ ( r ) πrdr = Z − λ I λ r ) πrdr = 8 π − λ/I λ (63)and I λ = π + 18 λ (64)lim λ →∞ λI λ = 8 . (65)The curve λI λ is shown on Fig. 2 and | u λ (1) − u λ (0) | in Fig. 4. Finally, u λ ( r ) = log(1 − λλ + 8 π r ) (66) | u λ (1) − u λ (0) | = 2 log(1 − λλ + 8 π ) . (67)We conclude that u λ ( r ) decreases smoothly and in the limit λ → ∞ , the solution blow-upover the entire boundary. 17 .3 Regular expansion of solution 15 for small λ We shall now study the small asymptotic expansion of solution 15 for small λ . Using aregular expansion, u ( x ) = u ( x ) + u ( x ) λ + u ( x ) λ + o ( λ ) , (68)we obtain using eq. 15 that u ( x ) = 0 and u is solution of − ∆ u = 1 | Ω | on Ω (69) ∂u ∂ n = − | ∂ Ω | on ∂ Ω . (70)For R = 1, u ( r ) = − r π . (71)with u (0) = 0. We conclude that u ( r ) ≤
0, Thus, u ( r ) = − r π λ + O ( λ ) . (72)The second order term u is solution of − ∆ u = − u | Ω | on Ω , (73)with u (0) = 0 and u ′ (0) = 0. For R = 1, u ( r ) = − r π . (74)Thus, u ( r ) = − r π λ − r π λ + O ( λ ) . (75)18 eferences [1] R. Araya, K.B. Eisenthal, R. Yuste, ”Dendritic spines linearize the summation ofexcitatory potentials”, Proc. Natl. Acad. Sci. USA
Proc. Natl. Acad. Sci. USA
Proc. Natl. Acad. Sci. USA , 104(30), 12347-52 (Jul., 2007).[4] R.A. Askey, Handbooks of special functions , Cambridge University Press; 1 edition(Feb. 15, 2001).[5] V. Barcilon, ”Ion Flow Through Narrow Membrane Channels: Part I”,
SIAM Journalon Applied Mathematics , (5), pp. 1391-1404 (Oct., 1992).[6] V. Barcilon, D.P. Chen, R.S. Eisenberg, ”Ion Flow through Narrow Membrane Chan-nels: Part II”, SIAM Journal on Applied Mathematics , (5), 1405-1425 (Oct.,1992).[7] A. Biess, E. Korkotian, D. Holcman, ”Diffusion in a dendritic spine: the role ofgeometry”, Phys. Rev. E. Stat. Nonlin. Soft Matter Phys (2 Pt 1), 021922. (Aug.2007).[8] B.L. Bloodgood, B.L. Sabatini, ”Neuronal activity regulates diffusion across the neckof dendritic spines”, Science , 310(5749), 866-9 (Nov., 2005).[9] W. Chen, R. Erban, S. Jonathan Chapman, ”From Brownian Dynamics to MarkovChain: An Ion Channel Example”, SIAM Journal of Applied Mathematics (1),208-235 (2014).[10] D.E. Goldman, ”Potential, impedance, and rectification in membranes”, J. Gen.Physiol. , (1), 37-60 (1943).[11] P. Graf, A. Nitzan, M.G. Kurnikova, R. Coalson, ”A Dynamic Lattice Monte CarloModel of Ion Transport in Inhomogeneous Dielectric Environments: Method andImplementation”, J. Phys. Chem. B , , 12324-12338 (2000).[12] B. Hille, ”Ionic Channels of Excitable Membranes 3rd edn”, Massachusetts: SinauerAssociates (2001).[13] D. Holcman, Z. Schuss, ”Time scale of diffusion in molecular and cellular biology”,
Journal of Physics A: Mathematical and Theoretical (17), 173001 (2014).[14] J. Jacobsen, K. Schmitt, ”The Liouville Bratu Gelfand Problem for radial operators, Journal of Differential Equations , pp.283-298 (2002).1915] J. Jacobsen, K. Schmitt, ”Radial solution of quasilinear elliptic differential equa-tions”,
Handbook of differential equations: Ordinary differential equations, 1st edn,Canada, Drabek and Fonda , chap. 4 (Sep., 2004).[16] D.D. Joseph T. S. Lundgren, ”Quasilinear Diricblet Problems Driven by PositiveSources”,
Archive for Rational Mechanics and Analysis , (4), pp 241-269 (1973).[17] C. Koch, T. Poggio, ”A Theoretical Analysis of Electrical Properties of Spines”, Proceedings of the Royal Society of London. Series B, Biological Sciences , (1213),pp. 455-477 (Jul. 22, 1983).[18] E. Korkotian, D. Holcman, M. Segal, ”Dynamic regulation of spine-dendrite couplingin cultured hippocampal neurons”,
Eur. J. Neurosci. (10), 2649-63 (Nov., 2004).[19] B. Nadler, Z. Schuss, A. Singer, R.S. Eisenberg, ”Ionic diffusion through confinedgeometries: from Langevin equations to partial differential equations”, J. Phys.:Condens. Matter , S2153-S2165 (2004).[20] Z. Schuss,”Theory and Applications of Stochastic Processes”, Springer series in Ap-plied Mathematical Sciences , Vol. 170, Springer Verlag, NY (2010).[21] Z. Schuss, B. Nadler, R.S. Eisenberg, ”Derivation of Poisson and Nernst-Planck equa-tions in a bath and channel from a molecular model”
Phys. Rev. E Stat. Nonlin. SoftMatter Phys. , 036116, (2001).[22] Z. Schuss, A. Singer, D. Holcman, ”The narrow escape problem for diffusion in cellularmicrodomains.” Proc. Natl. Acad. Sci. USA (41), pp.160981-161103 (2007).[23] I. Segev, W. Rall, ”Computational study of an excitable dendritic spine”
J Neuro-physiol , (2),499-523 (Aug., 1988).[24] A. Singer, Z. Schuss, ”Activation through a narrow opening”, SIAM J. Appl. Math. , (1), 98-108 (2007).[25] A. Singer, Z. Schuss, D. Holcman, ”NARROW ESCAPE, part III: Non-Smooth Do-mains and Riemann Surfaces”, J. Stat. Phys. ,(3), pp. 491-509 (2006).[26] K. Svoboda, D.W. Tank ,W. Denk, ”Direct measurement of coupling between den-dritic spines and shafts”,
Science , , 272(5262), 716-9 (1996).[27] A. Taflia, ”Diffusion of interacting particles in confined domains and applications tobiology”, Weizmann Institute, Phd dissertation, (2008).[28] V. H. Weston, ”On the asymptotic solution of a partial differential equation with anexponential nonlinearity”, SIAM Journal of Applied Mathematics , (6), (1978).2029] R. Yuste, W. Denk, ”Dendritic spines as basic functional units of neuronal integra-tion”, Nature , 375(6533), 682-4 (1995).[30] A. Zador ,C. Koch, T.H. Brown, ”Biophysical model of a Hebbian synapse”,
Proc.Natl. Acac. Sci. USA87