# On the Sodium Concentration Diffusion with Three-Dimensional Extracellular Stimulation

Authors: Luisa Consiglieri; Ana Rute Domingos
Journal: International Journal of Mathematics and Mathematical Sciences
 (2011)

---

## Abstract

We deal with the transmembrane sodium diffusion in a nerve. We study a mathematical model of a nerve fibre in response to an imposed extracellular stimulus. The presented model is constituted by a diffusion-drift vectorial equation in a bidomain, that is, two parabolic equations defined in each of the intra- and extra-regions. This system of partial differential equations can be understood as a reduced three-dimensional Poisson-Nernst-Planck model of the sodium concentration. The representation of the membrane includes a jump boundary condition describing the mechanisms involved in the excitation-contraction couple. Our first novelty comes from this general dynamical boundary condition. The second one is the three-dimensional behaviour of the extracellular stimulus. An analytical solution to the mathematical model is proposed depending on the morphology of the excitation.


---

## Body

## 1. Introduction

In the nervous system, there exists a cell transmembrane voltage due to the several types of ions on the opposite sites of the membrane [1, 2]. The ionic transmembrane flow is obtained by means of a given (mechanical, chemical, or electrical) signal. The action potential is generated on the membrane of the excitable cells. At the depolarisation stage, the inward sodium current appears if the voltage increases past a critical threshold, typically 15 mV higher than the resting value [3]. This runaway condition, whereby the positive feedback from the sodium current activates even more sodium channels, reveals the importance of the sodium ions among all ions presented in the axoplasm (the electrolytic fluid in the interior of the axon). Here, we deal with the ionic flow for the Na+ ions described by the ionic concentration (mol m−3). In order to understand the action potential and to offer predictions, the well-known Hodgkin-Huxley (HH for short) model plays an essential role for the quantitative understanding of the biological phenomena [4]. This work proposed that the action of potential in axon membranes can be analysed using cable theory. The authors proposed a system of four ordinary differential equations (ODEs) describing the current clamped experiments. Indeed, the previously unobserved dynamics in the HH model has a chaotic behaviour [5]. The field of computational neurophysiology has a long history containing extensive studies about the excitation of neural elements [6, 7]. A constructive discussion on the appropriate modelling of neural structures and their stimulation and blocking activities, by electrodes relatively remote from the target nerve cell, is provided in [6]. Rattay's book [7] illustrates whether the classical results for propagating action potentials, say the HH model for nonmyelinated fibres and the Frankenhauser-Huxley model for myelinated fibres, and subsequent analytical and numerical models may embody the phenomena and fit the electrophysical experiments. (The nerve cell or neuron is constituted by the soma (the cellular body), the dendrite, and one axon that connects the previous two. Some neurons have axons with an insulating layer, discontinuous, the so-called myelin sheath. These are the myelinated fibres. Neurons with naked axons, that is axons without myelin covering, are the so-called unmyelinated fibres [8].) Modified ODE systems [9–14] have extended the standard HH model and have been analysed through phase space methods (where equations are not explicitly solved). The control theory of the nonlinear systems exhibits chaotic behaviour of the version also known as the Fitzhugh-Nagumo (FN) model that consists of a second-order ODE dealing with the variation in time of the gating quantities and reinterprets the model developed by Hodgkin-Huxley [9]. Fitzhugh in [10] deals with a stable state and threshold phenomena as well as stable oscillations described by two variables of state, representing excitability and refractoriness, which are solutions of the so-called Bonhoeffer-van der Pol model. An extended FN system of ODE is numerically integrated in two different one-dimensional situations: free fibre and an externally stimulated clamped one [11]. A second-order differential equation of generalised FN type is solved by the least squares method, having as a solution the given single component (action potential) of the numerical solution [12]. Other variants of the HH model can be found in [13], based on geometric singular perturbation theory. Dynamics of spike initiation in other simplifications of the HH model, namely, the Morris-Lecar model, is exploited with phase plane and bifurcation analysis [14].Recently, using the Green function's method, analytic solutions for the cable equation response to the extracellular stimulus current have been found [15]. This alternative interpretation of the situation is the first step to understand the behaviour of the potential solution.Our concern is to understand the dynamics on the electrodiffusion of charged molecules (in particular, Na+ ions). We refer to [16, 17] for the 3D Poisson-Nernst-Planck (PNP) model analysed by finite element methods. Using a dynamic lattice Monte Carlo model [18], a description of the electrochemical processes is provided for the ion transport. After the reduction of the PNP model to a system of first-order ODE, in [19] the construction of singular orbits and the application of geometric singular perturbation theory provide information over permanently charged ions flow through an ion channel.Different mathematical models also study electrodiffusion. In [20], the authors look at a governing equation from the Maxwell equation in the quasistationary approximation when the electric potential jumps across the interface, and these jumps satisfy a dynamical condition (roughly speaking, in the form of a hyperbolic differential equation on the interface itself). In [21], the convergence of an electrochemical model is shown for a mixture of charged particles in a solution subject to prescribed electric potentials at two electrodes into a unique steady-state boundary value problem. An alternative approach, modelling the transmembrane potential in electrocardiology [22], considers a bidomain with a dynamic boundary jump condition, which is closely related to ours.The membrane dynamics used in this paper is based on the mathematical model started in the work [23]. Our model is derived from the Maxwell equations with current density defined by the Fick-Ohm law. Then, a drift term is included by the electrical contribution, which does not happen in the diffusion formulations obtained by mass and momentum conservation laws.Experimentally, an action potential is often generated by a rapid injection of current at a fixed point in the resting axon, which then spreads from the point of stimulation. However, the nerve cell is a three-dimensional structure. Even if a stimulus current pulse is arranged by the insertion of an electrode, a local current is developed. The membrane potential of the cell is not uniform at all points. The depolarisation spreading passively from an excited region of the membrane (near the insertion region of the electrode) to a neighbouring unexcited region occurs in three dimensions until uniformity is reached. This means that there is an interval of time where the flow has an angular dependence. The discrepancy between theory and experiments depends on the configuration of the experimental apparatus from which the propagated action potential was initiated and the strength of the current used to generate it [24]. Indeed, the discrepancy between the theoretical predicted and reported speeds of the propagated action potential is a consequence of neglecting the radial variation that occurs over small distances by comparison with axonal length. In [24], Fourier spectral methods are used to construct periodic solutions of the intracellular and extracellular potential for the Laplace equation.The goal of the present study is to determine how the profile of the activation affects the time parameter and the 3D domain of action potential initiating and, consequently, the propagation. This nonlinear model highlights the fact that an action potential is not generated instantaneously when the membrane potential crosses some preordained threshold [25]. We refer to [26] where the excitation response of an idealised infinite fibre is evaluated from the applied field of a unique point source electrode. Our main new contribution relies on the angular dependence of the concentration, since the stimulation can reliably propagate in a 3D form [25]. The description of the physiological phenomenon will be more realistic with 3D models. Moreover, we assume no azimuthal symmetry because it reflects the physical character of the biological phenomenon. In sum, we believe that these features can contribute to remove the actual discrepancy between the predicted and observed speed of the propagated action potential.Several biological constants are used throughout the paper. We keep them abstract so that the presented solution can be applied on different biological contexts. We illustrate their values with some examples (squid, cat, etc.).The paper is organised as follows. Section2 is devoted to state an initial and boundary value problem for the phenomenon under study. In Section 3 an analytical solution is obtained. Section 4 is a combination of results and some discussions about the presented model. In particular, a link to the now classical work of Hodgkin and Huxley (which can be applied to similar models) and how the executed technique can be useful to a particular example are discussed. Section 5 contains the conclusions. In the appendix, we briefly recall the theory for the confluent hypergeometric equation.


## 2. Statement of the Problem

The axon is a thin cellular extension, that may be short or long, responsible for transporting the information (electrical impulses) from the soma to the dendrite [1]. The axon can be described as a cylinder of length ℓ and radius h, surrounded by a membrane Γm of negligible thickness (the cylinder surface) and immersed in an extracellular medium (2.1)Ωe:={(x,y,z)∈R3:0<x<l,h2<y2+z2<r2},
for some 0<h<r (see Figure 1). Let Ω⊂ℝ3 be a neighbourhood of the membrane Γm defined as (2.2)Ω:={(x,y,z)∈R3:0<x<l,(h-ϵ)2<y2+z2<r2},
for some 0<ϵ<h. Thus Ωi:=Ω∖Ω¯e and Γm:=Ω∖(Ωi∪Ωe) denote the intracellular space and the membrane surface, respectively. Define the external boundary (2.3)Γ:=]0;l[×{(y,z)∈R2:y2+z2=r2}.
Let T>0. The problem under study is defined by the system of parabolic equations (for details see [23]) (2.4)∂ci∂t-D∇2ci+σiεci=0inΩi×]0,T[,∂ce∂t-D∇2ce+σeεce=0inΩe×]0,T[,
where D is the sodium diffusion coefficient (D=0.267×10-9 m2 s−1), [27]) and ε is the sodium permittivity (ε=6.4×10-10 F m−1). The instant of time t is measured in seconds, ∂/∂t denotes the time derivative, and ∇2 represents the Laplacian. Here, ce and ci denote the sodium concentration, respectively, in Ωe and in Ωi, measured in mol m−3. The electrical conductivity σs, with s∈{i,e}, is considered homogeneous and constant at each subdomain (extra- and intracellular domains), measured in S m−1 [28]: (2.5)σe=13inΩe,σi=53inΩi.
Other values for the diffusion coefficient and the electrical conductivity can be found in [29].Schematic representations of the axon and its extracellular medium (not in scale) in thex-axis direction from the input (x=0) to the output (x=ℓ) sites (this figure was produced by the function ParametricPlot3D of the software Mathematica 5.2 developed by Wolfram Research, and Microsoft Office Power Point).
(a)(b)We assume the following boundary conditions:(2.6)∇ce⋅nΓ=0on⁡Γ×]0,T[,(2.7)α∂∂t[c]+β[c]=D∇ce⋅non⁡Γm×]0,T[,
where nΓ and n denote the outward normal to Γ and Γm, respectively.  The insulating boundary condition in (2.6) represents the zero outflow. The interface condition in (2.7) governs the evolution of the discontinuity of the concentration, taking into account that the concentration jumps [c]:=ce-ci across the membrane satisfy a dynamical condition [20]. Due to many conducting channels the lipid axon membrane exhibits a capacitive/conducting behaviour. It separates internal and external conducting solutions. Such a gap between two conductors forms a significant electrical capacitor. In living cells, the ions lost via ionic channels by diffusion are returned by ionic pumps in order to overcome the electrochemical gradient [27]. We can distinguish three main states of the channel: open, closed, and inactive. The opening of those channels requires several gating events. As in [20, 23], the gating functions are assumed positive real constants.At the instantt=0, an external stimulus Φ, which can be electrical, mechanical, or chemical, is applied to depolarise the resting membrane. On the external region in the two-dimensional boundary {x=0} the following condition holds: (2.8)∇ce⋅n(0,y,z,0)=Φ(y,z),
where n=(1,0,0). The 1D propagation behaviour of an external stimulus has a peak in the axon initial segment induced by the real (or simulated) synaptic inputs [30, 31]. Our firing pattern plotted in Figure 2 corresponds to the 3D description. In the sequel, the 3D domain is considered defined in cylindrical coordinates, that is, the x-axis indicates the longitudinal distance along the length of the axon, ρ-axis the radial distance measured from the centre, and θ-axis the angular measure that performs the real three-dimensional feature of the axon behaviour. Then, the function Φ can be given by(2.9)Φ(ρ,θ)=(κ1ρ-1/2exp⁡(τρ)+κ2ρ1/2)cos⁡(θ2),
with κ1,κ2,τ∈ℝ, for every (ρ,θ)∈[h;r[×[-π;π[ (see Figure 2). The angular position θ=0 corresponds to the closeness to the external source for the interval of time that current redistribution is not concluded (see [7, page 154], e.g., for the monopolar electrode).Figure 2
Adimensional plot of the mappingΦ considered in (2.9), for (ρ,θ)∈[h;r[×[-π;π[. The plotted surface illustrates the qualitative behaviour as function of the radius ρ and the angle θ, at the input site (x=0) and the initial instant of time (t=0) (this figure  was produced by Mathematica 5.2, developed by Wolfram Research).The initial condition is assumed constituted by an averaged condition:(2.10)∫Ωici(⋅,0)=Ci,∫Ωece(⋅,0)=Ce,
with the values of Ci and Ce depending on the physiological data (e.g., in a cat motoneuron Ci=15 mol m−3 and Ce=150 mol m−3 (cf. [3]), and for values at the various axons we refer to [32, 33] and references therein).Finally, in order for the model to be accurate the Boltzmann principle must be satisfied. For this purpose, it will be sufficient for the validation of the effect of the concentration ratio through the membrane to be commensurable with the corresponding Nernst potential𝒩 in the resting state, that is, at instance t=0 and on Γm(2.11)(ceci)ρ=h,t=0=N:=exp⁡[FzNa+VdRTr],
where Vd represents the depolarisation voltage, F is the Faraday constant (F=9.649×104 C mol−1, [3]), R denotes the universal gas constant (R=8.314 J (mol·K)−1), and Tr denotes the reference absolute temperature, in Kelvin.


## 3. Analytical Solution

In this section we present an analytical solution for (2.4), for s=i,e. For the sake of simplicity in notations, we will omit the index s when we look only at the generic equation. When we also take into account the boundary and the initial conditions, we will denote with an index s the solution and all the constants involved, for s=i,e.In cylindrical coordinates, both equations in (2.4) read as follows: (3.1)∂c∂t-D(∂2c∂x2+∂2c∂ρ2+1ρ∂c∂ρ+1ρ2∂2c∂θ2)+σεc=0,
with (x,ρ,θ,t)∈]0,ℓ[×]h-ϵ,r[×]-π,π[×]0,T[. Indeed, the above system of parabolic equations is considered together with the set of additional restraints (2.6)–(2.11). The validation of the initial and the boundary conditions specifies the choice of the several constants involved in the characterisation of the solutions, namely, ξ, δ0,i, δ1,i, δ2,i,ι, and η (see Section 4). Thus, we get a solution (see Section 3.1, for details): (3.2)ce(x,ρ,θ,t)=1ξ(κ1ρ-1/2exp⁡[τρ]+κ2ρ1/2)exp⁡[ξx+(ξ2D-σeε)t]cos⁡(θ2),
with ξ given by (3.25), and (3.3)ci(x,ρ,θ,t):=(δ1,iρ-1/2exp⁡[ξx+(ξ+ι)ρ+(ξ2D-σeε)t]+(δ2,iρ-1/2exp⁡[(η+ι)ρ]+δ0,iρ1/2)exp⁡[ηx+(η2D-σiε)t])cos⁡(θ2),
where δ0,i, δ1,i, δ2,i,ι, and η are correlated such that (3.31)–(3.38) hold.


### 3.1. A Formal Derivation

Here we show an analytical solution for (2.4). First we study the differential equation (3.1); then we find an explicit solution for both the extra- and intracellular concentrations. Inspired by the Fourier method of separation of variables, a useful method for finding a solution of a partial differential equation with several variables, we look for a solution of the form (3.4)c(x,ρ,θ,t)=Y(x,ρ,t)Z(θ).
Consequently, from substituting this in (3.1) we obtain (3.5)-1D∂tYY+∂xxYY+∂ρρYY+1ρ∂ρYY+1ρ2Z′′Z-σDε=0,
where ∂ϰY:=∂Y/∂ϰ and ∂ϰϰ=∂2Y/∂ϰ2, for ϰ∈{t,x,ρ}. Then there exists λ∈ℝ such that (3.6)Z′′(θ)=λZ(θ),λ=[1D(∂tYY(x,ρ,t)+σε)-∂xxYY(x,ρ,t)-∂ρρYY(x,ρ,t)-1ρ∂ρYY(x,ρ,t)]ρ2.Remark 3.1.
The second-order ordinary differential equation forZ has the following solutions.
Ifλ<0, then there exist real constants d1, d2 such that
(3.7)Z(θ)=d1cos⁡(|λ|θ)+d2sin(|λ|θ);
if λ=0, then there exist real constants d1, d2 such that Z(θ)=d1θ+d2; if λ>0, then there exist real constants d1, d2 such that
(3.8)Z(θ)=d1exp⁡(λθ)+d2exp⁡(-λθ).
We observe that the validation of the phenomenal data yields the choice of λ<0 and bounded.By Remark3.1, λ is negative. Thus, the second-order ordinary differential equation for Z has the solution given by (3.7), where d1, d2 are real constants. Since we are interested in a nonnegative valued solution for (2.4), we assume that Y and Z are both nonnegative functions. Then, from (3.7) we get (3.9)Z(θ)=dcos⁡(|λ|θ),
with d nonnegative real constant. We will now analyse the second equation from (3.6), that is,(3.10)ρ2[∂xxY+∂ρρY-1D∂tY-σDεY]+ρ∂ρY+λY=0.
Inspired by the theory for the confluent hypergeometric equations (see the appendix), we look for a function of the form (3.11)Y(x,ρ,t)=ρ-|λ|u(x,ρ,t)+ρ|λ|v(x,t)
that satisfies (3.10), where (3.12)u(x,ρ,t)=∑j=0∞fj(x+ρ+pDt)jexp⁡[ax+bρ+qDt],
and a, b, p, and q are real constants to be chosen later. Solving (3.10) for ρ|λ|v, we get (3.13)∂xxv-1D∂tv-σDεv=0.
Thus, we choose (see [34]) (3.14)v(x,t)=v0exp⁡[ξx+(ξ2D-σε)t],
with ξ∈ℝ, and v0>0 in order to get positive solutions. Solving (3.10) for ρ-|λ|u, we obtain (3.15)∑j=0+∞[(a2+b2-q-σεD)fj⋅(x+ρ+pDt)j=02(j+1)(j+2)fj+2+(2(a+b)-p)(j+1)fj+1+(a2+b2-q-σεD)fj]⋅(x+ρ+pDt)j=0.


#### 3.1.1. Exterior Case

For simplicity, we will keep the unknown constantsλ, d, fj, p, a, b, q, v0, and ξ without the extracellular subscripts. In order to find the extracellular concentration ce, we use the Neumann boundary condition (2.6) in cylindrical coordinates. Thus, from (3.11)–(3.14) we get (3.16)∑j=0∞[(j+1)fj+1+(b-|λ|r)fj](x+r+pDt)j=-|λ|r2|λ|-1v0exp⁡[(ξ-a)x-br+(ξ2-q-σeεD)Dt].
Applying the Taylor formula in (x-x0), with x0=-(r+pDt), and denoting Uj=j!fj,  we have (3.17)Uj+1+(b-|λ|r)Uj=-|λ|r2|λ|-1v0(ξ-a)j×exp⁡[(a-ξ-b)r+((a-ξ)p+ξ2-q-σeεD)Dt].
Since there is no dependence in time, we get(3.18)(ξ-a)p+q=ξ2-σeεD.DenotingM:=b-|λ|/r and N:=-|λ|r2|λ|-1v0exp[(a-ξ-b)r], it follows that (3.19)Uj+1=-MUj+N(ξ-a)j,Uj+2=M2Uj+N(-M+ξ-a)(ξ-a)j.
From (3.15) it follows that(3.20)2Uj+2+(2(a+b)-p)Uj+1+(a2+b2-q-σeεD)Uj=0.
Thus, introducing (3.19) into (3.20) and using (3.18), after some calculations we obtain Uj=(ξ-a)jf0, with (3.21)f0=|λ|r2|λ|-2v0exp⁡[(a-ξ-b)r](2ξr+2|λ|-pr)(a-b)2+2(a-b)|λ|/r+2|λ|/r2+p(b+ξ-a-|λ|/r)-ξ2.
Thus, (3.22)ue(x,ρ,t)=f0exp⁡[ξx+(ξ-a+b)ρ+(ξ2-σeεD)Dt],
and consequently (3.23)Ye(x,ρ,t)=(f0ρ-|λ|exp⁡[(ξ-a+b)ρ]+v0ρ|λ|)⋅exp⁡[ξx+(ξ2D-σeε)t].Notice that it still remains to findλ, d, a, b, v0, and ξ. Using (2.8), (3.9), and (3.23), we obtain |λ|=1/2 (hence, λ=-1/4), d=κ1/(f0ξ), b-a=τ-ξ, and v0=f0κ2/κ1, concluding that (3.24)ce(x,ρ,θ,t)=1ξ(κ1ρ-1/2exp⁡[τρ]+κ2ρ1/2)exp⁡[ξx+(ξ2D-σeε)t]cos⁡θ2.Finally, the initial condition (2.10) yields(3.25)exp⁡[ξl]-1ξ(κ1∫hrρ1/2exp⁡[τρ]dρ+25κ2(r5/2-h5/2))=Ce4.
This means that ξ is well defined.


#### 3.1.2. Interior Case

We are interested in finding the intracellular concentrationci using the forms (3.12) and (3.14): (3.26)ci(x,ρ,θ,t)=(dρ-1/2u(x,ρ,t)+δ0,iρ1/2exp⁡[ηx+(η2D-σiε)t])cos⁡θ2,
with δ0,i=dv0=κ2/ξ and λ=-1/4. Here, we consider new unknown constants d, fj, p, a, b, q, δ0,i, and η (some of them relabelled to avoid the use of the intracellular subscripts).First, notice that (3.15) is verified under these new unknown constants, and again (3.15) implies (3.20) with σe replaced by σi, that is,(3.27)2Uj+2+(2(a+b)-p)Uj+1+(a2+b2-q-σiεD)Uj=0.
Next, we verify (2.7). Introducing (3.24) and (3.26) in (2.7), we obtain (3.28)∑j=0∞d(αpDfj+1(j+1)+(αqD+β)fj)(x+h+pDt)j=-δ0,ihexp⁡[-bh+(η-a)x+(η2D-σiε-qD)t](α(η2D-σiε)+β)+exp⁡[-bh+(ξ-a)x+(ξ2D-σe/ε-qD)t]ξ×(κ1(α(ξ2D-σeε)+β+D(12-hτ))exp⁡(τh)+κ2h(α(ξ2D-σeε)+β-D2)).
Arguing as in the exterior case (see Section 3.1.1), we can apply the Taylor formula in (x-x0) with x0=-(h+pDt) resulting in for each time level, (3.29)Uj+1=-MUj+N(ξ-a)j+P(η-a)j,Uj+2=M2Uj+N(-M+ξ-a)(ξ-a)j+P(-M+η-a)(η-a)j,
where M:=q/p+β/(αpD) and (3.30)P(t):=-δ0,ihdαpDexp⁡[(a-η-b)h+((a-η)p+η2-σiεD-q)Dt]×(α(η2D-σiε)+β),N(t):=1dαpDξexp⁡[(a-ξ-b)h+((a-ξ)p+ξ2-σeεD-q)Dt]×(κ1(α(ξ2D-σeε)+β+D(12-hτ))exp⁡(τh)+κ2h(α(ξ2D-σeε)+β-D2)).Introducing (3.29) into (3.27), we conclude that dUj=j!dfj=(ξ-a)jδ1,i+(η-a)jδ2,i, where (3.31)δ1,id=2NM-ξ-b+p/22M2-M(2(a+b)-p)+a2+b2-q-σi/(εD),(3.32)δ2,id=2PM-η-b+p/22M2-M(2(a+b)-p)+a2+b2-q-σi/(εD).
This implies the existence of the solution (3.3) if there is time independence in (3.31) and (3.32). This time independence can be given into (3.30) implying (3.33)(η-a)p+q=η2-σiεD,(3.34)(ξ-a)p+q=ξ2-σeεD.
These mean that the constants p and q are well defined, and P(t)≡P and N(t)≡N for all t.Notice that the constantsδ1,i, and δ2,i are well defined observing that the constants a, b, δ0,i and η will be known. For the sake of simplicity, if we take ι=b-a, then (3.26) leads us to (3.3). Thus, condition (2.10) reads (3.35)Ci4=exp⁡[ξl]-1ξδ1,i∫h-ϵhρ1/2exp⁡[(ξ+ι)ρ]dρ+exp⁡[ηl]-1η×(δ2,i∫h-ϵhρ1/2exp⁡[(η+ι)ρ]dρ+2δ0,i5(h5/2-(h-ϵ)5/2)).
This yields the choice for η. Finally, from (2.11) it follows that (3.36)1ξ(k1hexp⁡[τh]+k2h)exp⁡[ξx]=N(δ1,ihexp⁡[ξx+(ξ+ι)h]+δ2,ihexp⁡[ηx+(η+ι)h]+δ0,ihexp⁡[ηx]),
for all x∈]0,ℓ[. Hence, the algebraic system for ι and δ0,i is(3.37)1ξ(k1hexp⁡[τh]+k2h)=Nδ1,ihexp⁡[(ξ+ι)h],(3.38)δ0,i=-δ2,ihexp⁡[(η+ι)h],
concluding their existence.


## 3.1. A Formal Derivation

Here we show an analytical solution for (2.4). First we study the differential equation (3.1); then we find an explicit solution for both the extra- and intracellular concentrations. Inspired by the Fourier method of separation of variables, a useful method for finding a solution of a partial differential equation with several variables, we look for a solution of the form (3.4)c(x,ρ,θ,t)=Y(x,ρ,t)Z(θ).
Consequently, from substituting this in (3.1) we obtain (3.5)-1D∂tYY+∂xxYY+∂ρρYY+1ρ∂ρYY+1ρ2Z′′Z-σDε=0,
where ∂ϰY:=∂Y/∂ϰ and ∂ϰϰ=∂2Y/∂ϰ2, for ϰ∈{t,x,ρ}. Then there exists λ∈ℝ such that (3.6)Z′′(θ)=λZ(θ),λ=[1D(∂tYY(x,ρ,t)+σε)-∂xxYY(x,ρ,t)-∂ρρYY(x,ρ,t)-1ρ∂ρYY(x,ρ,t)]ρ2.Remark 3.1.
The second-order ordinary differential equation forZ has the following solutions.
Ifλ<0, then there exist real constants d1, d2 such that
(3.7)Z(θ)=d1cos⁡(|λ|θ)+d2sin(|λ|θ);
if λ=0, then there exist real constants d1, d2 such that Z(θ)=d1θ+d2; if λ>0, then there exist real constants d1, d2 such that
(3.8)Z(θ)=d1exp⁡(λθ)+d2exp⁡(-λθ).
We observe that the validation of the phenomenal data yields the choice of λ<0 and bounded.By Remark3.1, λ is negative. Thus, the second-order ordinary differential equation for Z has the solution given by (3.7), where d1, d2 are real constants. Since we are interested in a nonnegative valued solution for (2.4), we assume that Y and Z are both nonnegative functions. Then, from (3.7) we get (3.9)Z(θ)=dcos⁡(|λ|θ),
with d nonnegative real constant. We will now analyse the second equation from (3.6), that is,(3.10)ρ2[∂xxY+∂ρρY-1D∂tY-σDεY]+ρ∂ρY+λY=0.
Inspired by the theory for the confluent hypergeometric equations (see the appendix), we look for a function of the form (3.11)Y(x,ρ,t)=ρ-|λ|u(x,ρ,t)+ρ|λ|v(x,t)
that satisfies (3.10), where (3.12)u(x,ρ,t)=∑j=0∞fj(x+ρ+pDt)jexp⁡[ax+bρ+qDt],
and a, b, p, and q are real constants to be chosen later. Solving (3.10) for ρ|λ|v, we get (3.13)∂xxv-1D∂tv-σDεv=0.
Thus, we choose (see [34]) (3.14)v(x,t)=v0exp⁡[ξx+(ξ2D-σε)t],
with ξ∈ℝ, and v0>0 in order to get positive solutions. Solving (3.10) for ρ-|λ|u, we obtain (3.15)∑j=0+∞[(a2+b2-q-σεD)fj⋅(x+ρ+pDt)j=02(j+1)(j+2)fj+2+(2(a+b)-p)(j+1)fj+1+(a2+b2-q-σεD)fj]⋅(x+ρ+pDt)j=0.


### 3.1.1. Exterior Case

For simplicity, we will keep the unknown constantsλ, d, fj, p, a, b, q, v0, and ξ without the extracellular subscripts. In order to find the extracellular concentration ce, we use the Neumann boundary condition (2.6) in cylindrical coordinates. Thus, from (3.11)–(3.14) we get (3.16)∑j=0∞[(j+1)fj+1+(b-|λ|r)fj](x+r+pDt)j=-|λ|r2|λ|-1v0exp⁡[(ξ-a)x-br+(ξ2-q-σeεD)Dt].
Applying the Taylor formula in (x-x0), with x0=-(r+pDt), and denoting Uj=j!fj,  we have (3.17)Uj+1+(b-|λ|r)Uj=-|λ|r2|λ|-1v0(ξ-a)j×exp⁡[(a-ξ-b)r+((a-ξ)p+ξ2-q-σeεD)Dt].
Since there is no dependence in time, we get(3.18)(ξ-a)p+q=ξ2-σeεD.DenotingM:=b-|λ|/r and N:=-|λ|r2|λ|-1v0exp[(a-ξ-b)r], it follows that (3.19)Uj+1=-MUj+N(ξ-a)j,Uj+2=M2Uj+N(-M+ξ-a)(ξ-a)j.
From (3.15) it follows that(3.20)2Uj+2+(2(a+b)-p)Uj+1+(a2+b2-q-σeεD)Uj=0.
Thus, introducing (3.19) into (3.20) and using (3.18), after some calculations we obtain Uj=(ξ-a)jf0, with (3.21)f0=|λ|r2|λ|-2v0exp⁡[(a-ξ-b)r](2ξr+2|λ|-pr)(a-b)2+2(a-b)|λ|/r+2|λ|/r2+p(b+ξ-a-|λ|/r)-ξ2.
Thus, (3.22)ue(x,ρ,t)=f0exp⁡[ξx+(ξ-a+b)ρ+(ξ2-σeεD)Dt],
and consequently (3.23)Ye(x,ρ,t)=(f0ρ-|λ|exp⁡[(ξ-a+b)ρ]+v0ρ|λ|)⋅exp⁡[ξx+(ξ2D-σeε)t].Notice that it still remains to findλ, d, a, b, v0, and ξ. Using (2.8), (3.9), and (3.23), we obtain |λ|=1/2 (hence, λ=-1/4), d=κ1/(f0ξ), b-a=τ-ξ, and v0=f0κ2/κ1, concluding that (3.24)ce(x,ρ,θ,t)=1ξ(κ1ρ-1/2exp⁡[τρ]+κ2ρ1/2)exp⁡[ξx+(ξ2D-σeε)t]cos⁡θ2.Finally, the initial condition (2.10) yields(3.25)exp⁡[ξl]-1ξ(κ1∫hrρ1/2exp⁡[τρ]dρ+25κ2(r5/2-h5/2))=Ce4.
This means that ξ is well defined.


### 3.1.2. Interior Case

We are interested in finding the intracellular concentrationci using the forms (3.12) and (3.14): (3.26)ci(x,ρ,θ,t)=(dρ-1/2u(x,ρ,t)+δ0,iρ1/2exp⁡[ηx+(η2D-σiε)t])cos⁡θ2,
with δ0,i=dv0=κ2/ξ and λ=-1/4. Here, we consider new unknown constants d, fj, p, a, b, q, δ0,i, and η (some of them relabelled to avoid the use of the intracellular subscripts).First, notice that (3.15) is verified under these new unknown constants, and again (3.15) implies (3.20) with σe replaced by σi, that is,(3.27)2Uj+2+(2(a+b)-p)Uj+1+(a2+b2-q-σiεD)Uj=0.
Next, we verify (2.7). Introducing (3.24) and (3.26) in (2.7), we obtain (3.28)∑j=0∞d(αpDfj+1(j+1)+(αqD+β)fj)(x+h+pDt)j=-δ0,ihexp⁡[-bh+(η-a)x+(η2D-σiε-qD)t](α(η2D-σiε)+β)+exp⁡[-bh+(ξ-a)x+(ξ2D-σe/ε-qD)t]ξ×(κ1(α(ξ2D-σeε)+β+D(12-hτ))exp⁡(τh)+κ2h(α(ξ2D-σeε)+β-D2)).
Arguing as in the exterior case (see Section 3.1.1), we can apply the Taylor formula in (x-x0) with x0=-(h+pDt) resulting in for each time level, (3.29)Uj+1=-MUj+N(ξ-a)j+P(η-a)j,Uj+2=M2Uj+N(-M+ξ-a)(ξ-a)j+P(-M+η-a)(η-a)j,
where M:=q/p+β/(αpD) and (3.30)P(t):=-δ0,ihdαpDexp⁡[(a-η-b)h+((a-η)p+η2-σiεD-q)Dt]×(α(η2D-σiε)+β),N(t):=1dαpDξexp⁡[(a-ξ-b)h+((a-ξ)p+ξ2-σeεD-q)Dt]×(κ1(α(ξ2D-σeε)+β+D(12-hτ))exp⁡(τh)+κ2h(α(ξ2D-σeε)+β-D2)).Introducing (3.29) into (3.27), we conclude that dUj=j!dfj=(ξ-a)jδ1,i+(η-a)jδ2,i, where (3.31)δ1,id=2NM-ξ-b+p/22M2-M(2(a+b)-p)+a2+b2-q-σi/(εD),(3.32)δ2,id=2PM-η-b+p/22M2-M(2(a+b)-p)+a2+b2-q-σi/(εD).
This implies the existence of the solution (3.3) if there is time independence in (3.31) and (3.32). This time independence can be given into (3.30) implying (3.33)(η-a)p+q=η2-σiεD,(3.34)(ξ-a)p+q=ξ2-σeεD.
These mean that the constants p and q are well defined, and P(t)≡P and N(t)≡N for all t.Notice that the constantsδ1,i, and δ2,i are well defined observing that the constants a, b, δ0,i and η will be known. For the sake of simplicity, if we take ι=b-a, then (3.26) leads us to (3.3). Thus, condition (2.10) reads (3.35)Ci4=exp⁡[ξl]-1ξδ1,i∫h-ϵhρ1/2exp⁡[(ξ+ι)ρ]dρ+exp⁡[ηl]-1η×(δ2,i∫h-ϵhρ1/2exp⁡[(η+ι)ρ]dρ+2δ0,i5(h5/2-(h-ϵ)5/2)).
This yields the choice for η. Finally, from (2.11) it follows that (3.36)1ξ(k1hexp⁡[τh]+k2h)exp⁡[ξx]=N(δ1,ihexp⁡[ξx+(ξ+ι)h]+δ2,ihexp⁡[ηx+(η+ι)h]+δ0,ihexp⁡[ηx]),
for all x∈]0,ℓ[. Hence, the algebraic system for ι and δ0,i is(3.37)1ξ(k1hexp⁡[τh]+k2h)=Nδ1,ihexp⁡[(ξ+ι)h],(3.38)δ0,i=-δ2,ihexp⁡[(η+ι)h],
concluding their existence.


## 3.1.1. Exterior Case

For simplicity, we will keep the unknown constantsλ, d, fj, p, a, b, q, v0, and ξ without the extracellular subscripts. In order to find the extracellular concentration ce, we use the Neumann boundary condition (2.6) in cylindrical coordinates. Thus, from (3.11)–(3.14) we get (3.16)∑j=0∞[(j+1)fj+1+(b-|λ|r)fj](x+r+pDt)j=-|λ|r2|λ|-1v0exp⁡[(ξ-a)x-br+(ξ2-q-σeεD)Dt].
Applying the Taylor formula in (x-x0), with x0=-(r+pDt), and denoting Uj=j!fj,  we have (3.17)Uj+1+(b-|λ|r)Uj=-|λ|r2|λ|-1v0(ξ-a)j×exp⁡[(a-ξ-b)r+((a-ξ)p+ξ2-q-σeεD)Dt].
Since there is no dependence in time, we get(3.18)(ξ-a)p+q=ξ2-σeεD.DenotingM:=b-|λ|/r and N:=-|λ|r2|λ|-1v0exp[(a-ξ-b)r], it follows that (3.19)Uj+1=-MUj+N(ξ-a)j,Uj+2=M2Uj+N(-M+ξ-a)(ξ-a)j.
From (3.15) it follows that(3.20)2Uj+2+(2(a+b)-p)Uj+1+(a2+b2-q-σeεD)Uj=0.
Thus, introducing (3.19) into (3.20) and using (3.18), after some calculations we obtain Uj=(ξ-a)jf0, with (3.21)f0=|λ|r2|λ|-2v0exp⁡[(a-ξ-b)r](2ξr+2|λ|-pr)(a-b)2+2(a-b)|λ|/r+2|λ|/r2+p(b+ξ-a-|λ|/r)-ξ2.
Thus, (3.22)ue(x,ρ,t)=f0exp⁡[ξx+(ξ-a+b)ρ+(ξ2-σeεD)Dt],
and consequently (3.23)Ye(x,ρ,t)=(f0ρ-|λ|exp⁡[(ξ-a+b)ρ]+v0ρ|λ|)⋅exp⁡[ξx+(ξ2D-σeε)t].Notice that it still remains to findλ, d, a, b, v0, and ξ. Using (2.8), (3.9), and (3.23), we obtain |λ|=1/2 (hence, λ=-1/4), d=κ1/(f0ξ), b-a=τ-ξ, and v0=f0κ2/κ1, concluding that (3.24)ce(x,ρ,θ,t)=1ξ(κ1ρ-1/2exp⁡[τρ]+κ2ρ1/2)exp⁡[ξx+(ξ2D-σeε)t]cos⁡θ2.Finally, the initial condition (2.10) yields(3.25)exp⁡[ξl]-1ξ(κ1∫hrρ1/2exp⁡[τρ]dρ+25κ2(r5/2-h5/2))=Ce4.
This means that ξ is well defined.


## 3.1.2. Interior Case

We are interested in finding the intracellular concentrationci using the forms (3.12) and (3.14): (3.26)ci(x,ρ,θ,t)=(dρ-1/2u(x,ρ,t)+δ0,iρ1/2exp⁡[ηx+(η2D-σiε)t])cos⁡θ2,
with δ0,i=dv0=κ2/ξ and λ=-1/4. Here, we consider new unknown constants d, fj, p, a, b, q, δ0,i, and η (some of them relabelled to avoid the use of the intracellular subscripts).First, notice that (3.15) is verified under these new unknown constants, and again (3.15) implies (3.20) with σe replaced by σi, that is,(3.27)2Uj+2+(2(a+b)-p)Uj+1+(a2+b2-q-σiεD)Uj=0.
Next, we verify (2.7). Introducing (3.24) and (3.26) in (2.7), we obtain (3.28)∑j=0∞d(αpDfj+1(j+1)+(αqD+β)fj)(x+h+pDt)j=-δ0,ihexp⁡[-bh+(η-a)x+(η2D-σiε-qD)t](α(η2D-σiε)+β)+exp⁡[-bh+(ξ-a)x+(ξ2D-σe/ε-qD)t]ξ×(κ1(α(ξ2D-σeε)+β+D(12-hτ))exp⁡(τh)+κ2h(α(ξ2D-σeε)+β-D2)).
Arguing as in the exterior case (see Section 3.1.1), we can apply the Taylor formula in (x-x0) with x0=-(h+pDt) resulting in for each time level, (3.29)Uj+1=-MUj+N(ξ-a)j+P(η-a)j,Uj+2=M2Uj+N(-M+ξ-a)(ξ-a)j+P(-M+η-a)(η-a)j,
where M:=q/p+β/(αpD) and (3.30)P(t):=-δ0,ihdαpDexp⁡[(a-η-b)h+((a-η)p+η2-σiεD-q)Dt]×(α(η2D-σiε)+β),N(t):=1dαpDξexp⁡[(a-ξ-b)h+((a-ξ)p+ξ2-σeεD-q)Dt]×(κ1(α(ξ2D-σeε)+β+D(12-hτ))exp⁡(τh)+κ2h(α(ξ2D-σeε)+β-D2)).Introducing (3.29) into (3.27), we conclude that dUj=j!dfj=(ξ-a)jδ1,i+(η-a)jδ2,i, where (3.31)δ1,id=2NM-ξ-b+p/22M2-M(2(a+b)-p)+a2+b2-q-σi/(εD),(3.32)δ2,id=2PM-η-b+p/22M2-M(2(a+b)-p)+a2+b2-q-σi/(εD).
This implies the existence of the solution (3.3) if there is time independence in (3.31) and (3.32). This time independence can be given into (3.30) implying (3.33)(η-a)p+q=η2-σiεD,(3.34)(ξ-a)p+q=ξ2-σeεD.
These mean that the constants p and q are well defined, and P(t)≡P and N(t)≡N for all t.Notice that the constantsδ1,i, and δ2,i are well defined observing that the constants a, b, δ0,i and η will be known. For the sake of simplicity, if we take ι=b-a, then (3.26) leads us to (3.3). Thus, condition (2.10) reads (3.35)Ci4=exp⁡[ξl]-1ξδ1,i∫h-ϵhρ1/2exp⁡[(ξ+ι)ρ]dρ+exp⁡[ηl]-1η×(δ2,i∫h-ϵhρ1/2exp⁡[(η+ι)ρ]dρ+2δ0,i5(h5/2-(h-ϵ)5/2)).
This yields the choice for η. Finally, from (2.11) it follows that (3.36)1ξ(k1hexp⁡[τh]+k2h)exp⁡[ξx]=N(δ1,ihexp⁡[ξx+(ξ+ι)h]+δ2,ihexp⁡[ηx+(η+ι)h]+δ0,ihexp⁡[ηx]),
for all x∈]0,ℓ[. Hence, the algebraic system for ι and δ0,i is(3.37)1ξ(k1hexp⁡[τh]+k2h)=Nδ1,ihexp⁡[(ξ+ι)h],(3.38)δ0,i=-δ2,ihexp⁡[(η+ι)h],
concluding their existence.


## 4. Results and Discussion

Our continuum model (2.4)–(2.11) is developed to identify and analyse the diffusion phenomena, focusing on the average density distribution of species of charged particles and their description through unified partial differential equations (PDEs). We showed that these PDEs admit the existence of explicit solutions (3.2) and (3.3) for the intra- and extracellular sodium concentrations with well-determined constants, namely,(i)
κ1, κ2, and τ come from the profile of the extracellular stimulation (2.9),(ii)
ξ is given by (3.25), which is calculated from the initial condition (2.10) for the intracellular sodium concentration,(iii)
η is given by (3.35), which is calculated from the initial condition (2.10) for the extracellular sodium concentration,(iv)
ι and δ0,i are determined in (3.37) and (3.38), which come from the Boltzmann principle (2.11),(v)
δ1,i and δ2,i are determined in (3.31) and (3.32), which come from the jump condition (2.7).This result relates to the underlying mechanisms of electrodiffusion in sodium ions. However, we emphasise that it can be applied to any charged particles.The axon membrane acts as an interface between the intra- and extracellular concentrations of the Na+ ions where a dynamical jump condition is taken into account (see (2.7)). This avoids studying the ionic diffusion across the membrane. Indeed, the membrane is regarded such that the sodium concentration defined in the whole domain may suffer a finite jump when the membrane surface is crossed. The limit values of the sodium concentration may not, therefore, be the same when the membrane is approached from either the exterior or the interior. The concentration difference between the outside and inside of the membrane in (2.7) transports ions against their concentration gradients (from regions of high concentration to regions of low concentration). Therefore, this jump interface condition includes the continuous process of autotransformation of the membrane due to the movement of the proteins that transport ions. The first time derivative produces irreversible motion due to the exponential time-dependent factor. In order to determine the validity of the presence of α and β in (2.7), we pay attention to the HH model, as most of the models in electrophysiology are either its variants or its simplifications. The sodium conductance GNa, measured in S m−2, is regulated by voltage-dependent activation and inactivation variables usually called gating variables. The sodium conductance of the squid axon membrane is a product of three terms: a scale factor term gNa, a turning-on process m3, and an inactivating process h, where m and h denote the activation and the inactivation of the sodium current, respectively [3, 4, 7]: (4.1)GNa=gNam3h=120m3h.
Their dynamics are described by the ODE system: (4.2)dmdt=0.1(25-V)exp⁡[(25-V)/10]-1(1-m)-4exp⁡[V/18]m,dhdt=0.07exp⁡[V/20](1-h)-1exp⁡[(30-V)/10]+1h,
where V=Vm-VNa, with Vm and VNa representing the membrane potential and sodium equilibrium potential, in mV. This formulation expresses that the sodium conductance activation and inactivation are decoupled variables. A revised version of the HH model based on the molecular reaction sequence to account for both sodium and potassium conductance transients (Goldman-Hodgkin-Katz equation) was introduced by Clay [35]. Other kinetic cycles were proposed in [36] and references therein. Clearly, all these coincide with the HH model when one species is taken into account. Therefore, if the thickness of the membrane has a positive Lebesgue measure, the sodium conductance system can be modelled by expressing the kinetic model output as (4.3)I=Cm∂V∂t+(Vm-VNa)GNa,
with Cm denoting the membrane capacitance (F m−2). Since (4.2) has usually been solved under steady-state gating functions, which means that V≡V(x), then, passing to the limit as the thickness tends to zero, the steady-state conduction gives the characterisation (4.4)β~σm(membrane  conductivity),
considering the Poisson equation to compute the electric field from the charge present into the system (4.5)-∇⋅(ε∇V)=Fzc.
Using a time-dependent argument, α is correlated with the membrane permittivity. We refer to [37], a study of the electrical properties of the cell surface, known as the membrane.Also, our model fits the PNP model. Considering only one species, this electrodiffusion model is constituted by the Nernst-Planck equation to compute the ionic flux in an electrochemical gradient, in our notations,(4.6)∂c∂t+∇⋅D(∇c+FzkBTc∇V)=0,
and the adjoined Poisson equation (4.5), with kB representing the Boltzmann constant.Finally, we discuss how powerful is the technique executed here for the finding of explicit solutions to similar models. For instance, if we take the model from [38], with ci and ce given by (3.26) and (3.24), respectively, (4.7)∂ci∂t=αce-βci,
then equality (3.28) reads (4.8)∑j=0∞d(pDfj+1(j+1)+(qD+β(t))fj)(x+h+pDt)j=-α(t)δ0,ihexp⁡[-bh+(η-a)x+(η2D-σiε-qD)t]+exp⁡[-bh+(ξ-a)x+(ξ2D-σe/ε-qD)t]ξ×(κ1(α(t)+D(12-hτ))exp⁡(τh)+κ2h(α(t)-D2)),
observing that α≡α(t) and β≡β(t). Therefore, it follows that M(t):=q/p+β(t)/(pD),(4.9)P(t):=-α(t)δ0,ihdpDexp⁡[(a-η-b)h+((a-η)p+η2-σiεD-q)Dt],N(t):=1dpDξexp⁡[(a-ξ-b)h+((a-ξ)p+ξ2-σeεD-q)Dt]×(κ1(α(t)+D(12-hτ))exp⁡(τh)+κ2h(α(t)-D2)).
As in Section 3.1, expressions (3.31) and (3.32), with these new expressions, imply the existence of the solution (3.3) if there is independence in time. Thus, new relations can be obtained between the unknown constants of the executed technique and the physiological data, observing that the gating functions have exponential forms. Notice that the extracellular concentration keeps its definition, and this new jump condition (4.7) infers new constants into the definition for the intracellular concentration.The mobility of sodium channels, the mechanisms by which they are distributed over neurons, and the maintenance of this distribution are important issues for the physiology of excitable cells. Recent progress in determining 3D structures of biomolecules such as ion channels greatly facilitates the diffusion continuum description; see, for instance, [39]. The theory behind electrodynamics turns into the electrobiology by studying the diffusing channels which are simply reflected by the boundary.


## 5. Conclusions

The present work addresses the 3D effects of an initial stimulus applied on a section of the axon instead of over a point source. The obtained analytical solution demonstrates that the angularly distributed nature plays a crucial role in determining the action potential initiation. The sinusoidal shape over the angular structure is not documented in the literature since the experimental studies are devoted to the radial and longitudinal behaviours [40]. In [28], a 3D model predicts changes in the effects of the activation and inactivation gates of the sodium channel and consequently in the response on the action potential from the applied point source (in particular, the electrode position) relative to the geometry of the neuron.Here, the electrical conduction on the initial segment has the same pattern whether or not the axon is ensheathed in myelin [41]. Since the myelin sheath can be considered as a perfect insulator due to the core-conductor theory, the ionic flux only crosses the fibre at the nodes of Ranvier situated at the membrane of a myelinated axon. Therefore, our model can consider either unmyelinated or myelinated fibres taking the subunit Ω at each node compartment with ℓ representing the nodal gap width, that is, (2.7) describes the membrane dynamics in discrete space intervals. Even when an unmyelinated fibre is partly covered by Schwann cells, the axon membrane is separated from these cells by a space that is connected with the extracellular space [7, page 34].Our main result is the existence of explicit solutions for the intra- and extracellular sodium concentrations. From Section4, we can conclude that although the gating parameters were assumed constants our technique is sufficiently powerful to be applied to time-dependent parameters. Consequently, our model allows more complex solutions in accordance with the same physiological data.The existence of an analytical solution to our model constitutes the basis of ongoing numerical analysis, in order to understand the accurate evaluation during diffusion and to confirm the theoretical findings.



---

*Source: 862813-2011-07-21.xml*