On the Sodium Concentration Diffusion with Three-Dimensional Extracellular Stimulation

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.

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⁻³). 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. Section 2 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 ℎ, surrounded by a membrane Γ_𝑚 of negligible thickness (the cylinder surface) and immersed in an extracellular medium Ω_𝑒:=𝑥,𝑦,𝑧∈ℝ³:0<𝑥<ℓ,ℎ²<𝑦²+𝑧²<𝑟², (2.1) for some 0<ℎ<𝑟 (see Figure 1). Let Ω⊂ℝ³ be a neighbourhood of the membrane Γ_𝑚 defined as Ω:=𝑥,𝑦,𝑧∈ℝ³:0<𝑥<ℓ,(ℎ−𝜖)²<𝑦²+𝑧²<𝑟², (2.2) for some 0<𝜖<ℎ. Thus Ω_𝑖:=Ω⧵Ω_𝑒 and Γ_𝑚:=Ω⧵(Ω_𝑖∪Ω_𝑒) denote the intracellular space and the membrane surface, respectively. Define the external boundary Γ:=]0;ℓ[×𝑦,𝑧∈ℝ²:𝑦²+𝑧²=𝑟². (2.3) Let 𝑇>0. The problem under study is defined by the system of parabolic equations (for details see [23]) 𝜕𝑐_𝑖𝜕𝑡−𝐷∇²𝑐_𝑖+𝜎_𝑖𝜀𝑐_𝑖=0inΩ_𝑖×]0,𝑇[,𝜕𝑐_𝑒𝜕𝑡−𝐷∇²𝑐_𝑒+𝜎_𝑒𝜀𝑐_𝑒=0inΩ_𝑒×]0,𝑇[, (2.4) where 𝐷 is the sodium diffusion coefficient (𝐷=0.267×10⁻⁹ m² s⁻¹), [27]) and 𝜀 is the sodium permittivity (𝜀=6.4×10⁻¹⁰ F m⁻¹). The instant of time 𝑡 is measured in seconds, 𝜕/𝜕𝑡 denotes the time derivative, and ∇² represents the Laplacian. Here, 𝑐_𝑒 and 𝑐_𝑖 denote the sodium concentration, respectively, in Ω_𝑒 and in Ω_𝑖, measured in mol m⁻³. The electrical conductivity 𝜎_𝑠, with 𝑠∈{𝑖,𝑒}, is considered homogeneous and constant at each subdomain (extra- and intracellular domains), measured in S m⁻¹ [28]: 𝜎_𝑒=13inΩ_𝑒,𝜎_𝑖=53inΩ_𝑖. (2.5) Other values for the diffusion coefficient and the electrical conductivity can be found in [29].

(a)

(b)

Figure 1 

Schematic representations of the axon and its extracellular medium (not in scale) in the 𝑥-axis direction from the input (𝑥=0) to the output (𝑥=ℓ) sites (this figure was produced by the function ParametricPlot3D of the software Mathematica 5.2 developed by Wolfram Research, and Microsoft Office Power Point).

We assume the following boundary conditions:∇𝑐_𝑒⋅𝐧_Γ][,𝛼𝜕=0onΓ×0,𝑇[𝑐][𝑐]𝜕𝑡+𝛽=𝐷∇𝑐_𝑒⋅𝐧onΓ_𝑚×]0,𝑇[, (2.6)(2.7) where 𝐧_Γ and 𝐧 denote the outward normal to Γ and Γ_𝑚, 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 [𝑐]:=𝑐_𝑒−𝑐_𝑖 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 instant 𝑡=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 {𝑥=0} the following condition holds: ∇𝑐_𝑒⋅𝐧0,𝑦,𝑧,0=Φ𝑦,𝑧, (2.8) where 𝐧=(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 𝑥-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Φ𝜌,𝜃=𝜅₁𝜌^−1/2exp𝜏𝜌+𝜅₂𝜌^1/2𝜃cos2, (2.9) with 𝜅₁,𝜅₂,𝜏∈ℝ, for every (𝜌,𝜃)∈[ℎ;𝑟[×[−𝜋;𝜋[ (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 (𝜌,𝜃)∈[ℎ;𝑟[×[−𝜋;𝜋[. The plotted surface illustrates the qualitative behaviour as function of the radius 𝜌 and the angle 𝜃, at the input site (𝑥=0) and the initial instant of time (𝑡=0) (this figure was produced by Mathematica 5.2, developed by Wolfram Research).

The initial condition is assumed constituted by an averaged condition: _Ω𝑖𝑐_𝑖(⋅,0)=𝐶_𝑖,_Ω𝑒𝑐_𝑒(⋅,0)=𝐶_𝑒, (2.10) with the values of 𝐶_𝑖 and 𝐶_𝑒 depending on the physiological data (e.g., in a cat motoneuron 𝐶_𝑖=15 mol m⁻³ and 𝐶_𝑒=150 mol m⁻³ (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 𝑡=0 and on Γ_𝑚𝑐_𝑒𝑐_{𝑖𝜌=ℎ,𝑡=0}=𝒩:=exp𝐹𝑧_Na+𝑉_𝑑𝑅𝑇_𝑟, (2.11) where 𝑉_𝑑 represents the depolarisation voltage, 𝐹 is the Faraday constant (𝐹=9.649×10⁴ C mol⁻¹, [3]), 𝑅 denotes the universal gas constant (𝑅=8.314 J (mol·K)⁻¹), and 𝑇_𝑟 denotes the reference absolute temperature, in Kelvin.

3. Analytical Solution

In this section we present an analytical solution for (2.4), for 𝑠=𝑖,𝑒. For the sake of simplicity in notations, we will omit the index 𝑠 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 𝑠 the solution and all the constants involved, for 𝑠=𝑖,𝑒.

In cylindrical coordinates, both equations in (2.4) read as follows: 𝜕𝑐𝜕𝜕𝑡−𝐷²𝑐𝜕𝑥²+𝜕²𝑐𝜕𝜌²+1𝜌𝜕𝑐+1𝜕𝜌𝜌²𝜕²𝑐𝜕𝜃²+𝜎𝜀𝑐=0, (3.1) with (𝑥,𝜌,𝜃,𝑡)∈]0,ℓ[×]ℎ−𝜖,𝑟[×]−𝜋,𝜋[×]0,𝑇[. 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,𝑖, 𝛿_1,𝑖, 𝛿_2,𝑖,𝜄, and 𝜂 (see Section 4). Thus, we get a solution (see Section 3.1, for details): 𝑐_𝑒=1𝑥,𝜌,𝜃,𝑡𝜉𝜅₁𝜌^−1/2exp𝜏𝜌+𝜅₂𝜌^1/2exp𝜉𝑥+𝜉²𝜎𝐷−^𝑒𝜀𝑡𝜃cos2, (3.2) with 𝜉 given by (3.25), and 𝑐_𝑖𝑥,𝜌,𝜃,𝑡:=𝛿_1,𝑖𝜌^−1/2exp𝜉𝑥+(𝜉+𝜄)𝜌+𝜉²𝜎𝐷−^𝑒𝜀𝑡+𝛿_2,𝑖𝜌^−1/2exp𝜂+𝜄𝜌+𝛿_0,𝑖𝜌^1/2exp𝜂𝑥+𝜂²𝜎𝐷−^𝑖𝜀𝑡𝜃cos2, (3.3) where 𝛿_0,𝑖, 𝛿_1,𝑖, 𝛿_2,𝑖,𝜄, 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) Consequently, from substituting this in (3.1) we obtain −1𝐷𝜕_𝑡𝑌𝑌+𝜕_𝑥𝑥𝑌𝑌+𝜕_𝜌𝜌𝑌𝑌+1𝜌𝜕_𝜌𝑌𝑌+1𝜌²𝑍^′′𝑍−𝜎𝐷𝜀=0, (3.5) where 𝜕_𝜘𝑌:=𝜕𝑌/𝜕𝜘 and 𝜕_𝜘𝜘=𝜕²𝑌/𝜕𝜘², for 𝜘∈{𝑡,𝑥,𝜌}. Then there exists 𝜆∈ℝ such that 𝑍^′′(𝜃)(𝜃),1=𝜆𝑍𝜆=𝐷𝜕_𝑡𝑌𝑌+𝜎𝑥,𝜌,𝑡𝜀−𝜕_𝑥𝑥𝑌𝑌−𝜕𝑥,𝜌,𝑡^𝜌𝜌𝑌𝑌−1𝑥,𝜌,𝑡𝜌𝜕_𝜌𝑌𝑌𝑥,𝜌,𝑡𝜌². (3.6)

Remark 3.1. The second-order ordinary differential equation for Z has the following solutions.
If 𝜆<0, then there exist real constants 𝑑₁, 𝑑₂ such that 𝑍(𝜃)=𝑑₁cos|𝜆|𝜃+𝑑₂sin|𝜆|𝜃; (3.7) if 𝜆=0, then there exist real constants 𝑑₁, 𝑑₂ such that 𝑍(𝜃)=𝑑₁𝜃+𝑑₂; if 𝜆>0, then there exist real constants 𝑑₁, 𝑑₂ such that 𝑍(𝜃)=𝑑₁exp√𝜆𝜃+𝑑₂exp−√𝜆𝜃. (3.8) We observe that the validation of the phenomenal data yields the choice of 𝜆<0 and bounded.

By Remark 3.1, 𝜆 is negative. Thus, the second-order ordinary differential equation for 𝑍 has the solution given by (3.7), where 𝑑₁, 𝑑₂ are real constants. Since we are interested in a nonnegative valued solution for (2.4), we assume that 𝑌 and 𝑍 are both nonnegative functions. Then, from (3.7) we get 𝑍(𝜃)=𝑑cos|𝜆|𝜃, (3.9) with 𝑑 nonnegative real constant. We will now analyse the second equation from (3.6), that is,𝜌²𝜕_𝑥𝑥𝑌+𝜕_𝜌𝜌1𝑌−𝐷𝜕_𝑡𝜎𝑌−𝑌𝐷𝜀+𝜌𝜕_𝜌𝑌+𝜆𝑌=0. (3.10) Inspired by the theory for the confluent hypergeometric equations (see the appendix), we look for a function of the form 𝑌𝑥,𝜌,𝑡=𝜌^−|𝜆|𝑢𝑥,𝜌,𝑡+𝜌^|𝜆|𝑣(𝑥,𝑡) (3.11) that satisfies (3.10), where 𝑢𝑥,𝜌,𝑡=^∞_𝑗=0𝑓_𝑗𝑥+𝜌+𝑝𝐷𝑡^𝑗exp𝑎𝑥+𝑏𝜌+𝑞𝐷𝑡, (3.12) and 𝑎, 𝑏, 𝑝, and 𝑞 are real constants to be chosen later. Solving (3.10) for 𝜌^|𝜆|𝑣, we get 𝜕_𝑥𝑥1𝑣−𝐷𝜕_𝑡𝜎𝑣−𝐷𝜀𝑣=0. (3.13) Thus, we choose (see [34]) 𝑣(𝑥,𝑡)=𝑣₀exp𝜉𝑥+𝜉²𝜎𝐷−𝜀𝑡, (3.14) with 𝜉∈ℝ, and 𝑣₀>0 in order to get positive solutions. Solving (3.10) for 𝜌^−|𝜆|𝑢, we obtain ^+∞_𝑗=02𝑗+1𝑗+2𝑓_𝑗+2+2(𝑎+𝑏)−𝑝𝑗+1𝑓_𝑗+1+𝑎²+𝑏²𝜎−𝑞−𝑓𝜀𝐷^𝑗⋅𝑥+𝜌+𝑝𝐷𝑡^𝑗=0. (3.15)

3.1.1. Exterior Case

For simplicity, we will keep the unknown constants 𝜆, 𝑑, 𝑓_𝑗, 𝑝, 𝑎, 𝑏, 𝑞, 𝑣₀, and 𝜉 without the extracellular subscripts. In order to find the extracellular concentration 𝑐_𝑒, we use the Neumann boundary condition (2.6) in cylindrical coordinates. Thus, from (3.11)–(3.14) we get 𝑗+1𝑓^∞_𝑗=0𝑗+1+𝑏−|𝜆|𝑟𝑓_𝑗𝑥+𝑟+𝑝𝐷𝑡^𝑗=−|𝜆|𝑟^2|𝜆|−1𝑣₀exp(𝜉−𝑎)𝑥−𝑏𝑟+𝜉²𝜎−𝑞−^𝑒.𝜀𝐷𝐷𝑡 (3.16) Applying the Taylor formula in (𝑥−𝑥₀), with 𝑥₀=−(𝑟+𝑝𝐷𝑡), and denoting 𝑈_𝑗=𝑗!𝑓_𝑗, we have 𝑈_𝑗+1+𝑏−|𝜆|𝑟𝑈_𝑗=−|𝜆|𝑟^2|𝜆|−1𝑣₀(𝜉−𝑎)^𝑗×exp(𝑎−𝜉−𝑏)𝑟+(𝑎−𝜉)𝑝+𝜉²𝜎−𝑞−^𝑒.𝜀𝐷𝐷𝑡 (3.17) Since there is no dependence in time, we get(𝜉−𝑎)𝑝+𝑞=𝜉²−𝜎_𝑒.𝜀𝐷 (3.18)

Denoting 𝑀:=𝑏−|𝜆|/𝑟 and 𝑁:=−|𝜆|𝑟^2|𝜆|−1𝑣₀exp[(𝑎−𝜉−𝑏)𝑟], it follows that 𝑈_𝑗+1=−𝑀𝑈_𝑗+𝑁(𝜉−𝑎)^𝑗,𝑈_𝑗+2=𝑀²𝑈_𝑗+𝑁(−𝑀+𝜉−𝑎)(𝜉−𝑎)^𝑗. (3.19) From (3.15) it follows that2𝑈_𝑗+2+2(𝑎+𝑏)−𝑝𝑈_𝑗+1+𝑎²+𝑏²𝜎−𝑞−^𝑒𝑈𝜀𝐷^𝑗=0. (3.20) Thus, introducing (3.19) into (3.20) and using (3.18), after some calculations we obtain 𝑈_𝑗=(𝜉−𝑎)^𝑗𝑓₀, with 𝑓₀=|𝜆|𝑟^2|𝜆|−2𝑣₀exp(𝑎−𝜉−𝑏)𝑟2𝜉𝑟+2|𝜆|−𝑝𝑟(𝑎−𝑏)²+2(𝑎−𝑏)|𝜆|/𝑟+2|𝜆|/𝑟²+𝑝𝑏+𝜉−𝑎−|𝜆|/𝑟−𝜉². (3.21) Thus, 𝑢_𝑒𝑥,𝜌,𝑡=𝑓₀exp𝜉𝑥+(𝜉−𝑎+𝑏)𝜌+𝜉²−𝜎_𝑒,𝜀𝐷𝐷𝑡 (3.22) and consequently 𝑌_𝑒𝑥,𝜌,𝑡=𝑓₀𝜌^−|𝜆|exp(𝜉−𝑎+𝑏)𝜌+𝑣₀𝜌^|𝜆|⋅exp𝜉𝑥+𝜉²𝜎𝐷−^𝑒𝜀𝑡. (3.23)

Notice that it still remains to find 𝜆, 𝑑, 𝑎, 𝑏, 𝑣₀, and 𝜉. Using (2.8), (3.9), and (3.23), we obtain |𝜆|=1/2 (hence, 𝜆=−1/4), 𝑑=𝜅₁/(𝑓₀𝜉), 𝑏−𝑎=𝜏−𝜉, and 𝑣₀=𝑓₀𝜅₂/𝜅₁, concluding that 𝑐_𝑒=1𝑥,𝜌,𝜃,𝑡𝜉𝜅₁𝜌^−1/2exp𝜏𝜌+𝜅₂𝜌^1/2exp𝜉𝑥+𝜉²𝜎𝐷−^𝑒𝜀𝑡𝜃cos2. (3.24)

Finally, the initial condition (2.10) yieldsexp𝜉ℓ−1𝜉𝜅₁^𝑟_ℎ𝜌^1/22exp𝜏𝜌𝑑𝜌+5𝜅₂𝑟^5/2−ℎ^5/2=𝐶_𝑒4. (3.25) This means that 𝜉 is well defined.

3.1.2. Interior Case

We are interested in finding the intracellular concentration 𝑐_𝑖 using the forms (3.12) and (3.14): 𝑐_𝑖𝑥,𝜌,𝜃,𝑡=𝑑𝜌^−1/2𝑢𝑥,𝜌,𝑡+𝛿_0,𝑖𝜌^1/2exp𝜂𝑥+𝜂²𝜎𝐷−^𝑖𝜀𝑡𝜃cos2, (3.26) with 𝛿_0,𝑖=𝑑𝑣₀=𝜅₂/𝜉 and 𝜆=−1/4. Here, we consider new unknown constants 𝑑, 𝑓_𝑗, 𝑝, 𝑎, 𝑏, 𝑞, 𝛿_0,𝑖, 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 𝜎_𝑒 replaced by 𝜎_𝑖, that is,2𝑈_𝑗+2+2(𝑎+𝑏)−𝑝𝑈_𝑗+1+𝑎²+𝑏²𝜎−𝑞−^𝑖𝑈𝜀𝐷^𝑗=0. (3.27) Next, we verify (2.7). Introducing (3.24) and (3.26) in (2.7), we obtain ^∞_𝑗=0𝑑𝛼𝑝𝐷𝑓_𝑗+1𝑗+1+𝛼𝑞𝐷+𝛽𝑓_𝑗𝑥+ℎ+𝑝𝐷𝑡^𝑗=−𝛿_0,𝑖ℎexp−𝑏ℎ+𝜂−𝑎𝑥+𝜂²𝜎𝐷−^𝑖𝜀−𝑞𝐷𝑡𝛼𝜂²𝜎𝐷−^𝑖𝜀+()𝜉+𝛽exp−𝑏ℎ+𝜉−𝑎𝑥+²𝐷−𝜎_𝑒/𝜀−𝑞𝐷𝑡𝜉×𝜅₁𝛼𝜉²𝜎𝐷−^𝑒𝜀1+𝛽+𝐷2−ℎ𝜏exp(𝜏ℎ)+𝜅₂ℎ𝛼𝜉²𝜎𝐷−^𝑒𝜀𝐷+𝛽−2. (3.28) Arguing as in the exterior case (see Section 3.1.1), we can apply the Taylor formula in (𝑥−𝑥₀) with 𝑥₀=−(ℎ+𝑝𝐷𝑡) resulting in for each time level, 𝑈_𝑗+1=−𝑀𝑈_𝑗+𝑁(𝜉−𝑎)^𝑗+𝑃𝜂−𝑎^𝑗,𝑈_𝑗+2=𝑀²𝑈_𝑗+𝑁(−𝑀+𝜉−𝑎)(𝜉−𝑎)^𝑗+𝑃−𝑀+𝜂−𝑎𝜂−𝑎^𝑗, (3.29) where 𝑀∶=𝑞/𝑝+𝛽/(𝛼𝑝𝐷) and 𝑃(𝑡)𝛿:=−^0,𝑖ℎ𝑑𝛼𝑝𝐷exp𝑎−𝜂−𝑏ℎ+𝑎−𝜂𝑝+𝜂²−𝜎_𝑖×𝛼𝜂𝜀𝐷−𝑞𝐷𝑡₂𝜎𝐷−^𝑖𝜀,𝑁(𝑡)1+𝛽:=)()𝑑𝛼𝑝𝐷𝜉exp(𝑎−𝜉−𝑏ℎ+𝑎−𝜉𝑝+𝜉²−𝜎_𝑒×𝜅𝜀𝐷−𝑞𝐷𝑡₁𝛼𝜉²𝜎𝐷−^𝑒𝜀1+𝛽+𝐷2−ℎ𝜏exp(𝜏ℎ)+𝜅₂ℎ𝛼𝜉²𝜎𝐷−^𝑒𝜀𝐷+𝛽−2. (3.30)

Introducing (3.29) into (3.27), we conclude that 𝑑𝑈_𝑗=𝑗!𝑑𝑓_𝑗=(𝜉−𝑎)^𝑗𝛿_1,𝑖+(𝜂−𝑎)^𝑗𝛿_2,𝑖, where 𝛿_1,𝑖𝑑=2𝑁𝑀−𝜉−𝑏+𝑝/22𝑀²−𝑀2(𝑎+𝑏)−𝑝+𝑎²+𝑏²−𝑞−𝜎_𝑖/(),𝛿𝜀𝐷_2,𝑖𝑑=2𝑃𝑀−𝜂−𝑏+𝑝/22𝑀²−𝑀2(𝑎+𝑏)−𝑝+𝑎²+𝑏²−𝑞−𝜎_𝑖/().𝜀𝐷 (3.31)(3.32) 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)(3.34) These mean that the constants 𝑝 and 𝑞 are well defined, and 𝑃(𝑡)≡𝑃 and 𝑁(𝑡)≡𝑁 for all 𝑡.

Notice that the constants 𝛿_1,𝑖, and 𝛿_2,𝑖 are well defined observing that the constants 𝑎, 𝑏, 𝛿_0,𝑖 and 𝜂 will be known. For the sake of simplicity, if we take 𝜄=𝑏−𝑎, then (3.26) leads us to (3.3). Thus, condition (2.10) reads 𝐶_𝑖4=exp𝜉ℓ−1𝜉𝛿_1,𝑖^ℎ_ℎ−𝜖𝜌^1/2exp(𝜉+𝜄)𝜌𝑑𝜌+exp𝜂ℓ−1𝜂×𝛿_2,𝑖^ℎ_ℎ−𝜖𝜌^1/2exp𝜂+𝜄𝜌𝑑𝜌+2𝛿_0,𝑖5ℎ^5/2−(ℎ−𝜖)^5/2. (3.35) This yields the choice for 𝜂. Finally, from (2.11) it follows that 1𝜉𝑘₁√ℎexp[𝜏ℎ]+𝑘₂√ℎ𝛿exp𝜉𝑥=𝒩^1,𝑖√ℎ()ℎ+𝛿exp𝜉𝑥+𝜉+𝜄^2,𝑖√ℎexp𝜂𝑥+𝜂+𝜄ℎ+𝛿_0,𝑖√ℎexp𝜂𝑥, (3.36) for all 𝑥∈]0,ℓ[. Hence, the algebraic system for 𝜄 and 𝛿_0,𝑖 is1𝜉𝑘₁√ℎexp[𝜏ℎ]+𝑘₂√ℎ𝛿=𝒩^1,𝑖√ℎ)ℎ,𝛿exp(𝜉+𝜄_0,𝑖𝛿=−^2,𝑖ℎexp𝜂+𝜄ℎ, (3.37)(3.38) 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)𝜅₁, 𝜅₂, 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,𝑖 are determined in (3.37) and (3.38), which come from the Boltzmann principle (2.11),(v)𝛿_1,𝑖 and 𝛿_2,𝑖 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 𝐺_Na, measured in S m⁻², 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 𝑔_Na, a turning-on process 𝑚³, and an inactivating process ℎ, where 𝑚 and ℎ denote the activation and the inactivation of the sodium current, respectively [3, 4, 7]: 𝐺_Na=𝑔_Na𝑚³ℎ=120𝑚³ℎ. (4.1) Their dynamics are described by the ODE system: 𝑑𝑚=()𝑑𝑡0.125−𝑉)]()−4exp[(25−𝑉/10−11−𝑚exp[𝑉/18]𝑚,𝑑ℎ=𝑑𝑡0.07[]()−1exp𝑉/201−ℎexp[(30−𝑉)/10]+1ℎ, (4.2) where 𝑉=𝑉_𝑚−𝑉_Na, with 𝑉_𝑚 and 𝑉_Na 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 𝐼=𝐶_𝑚𝜕𝑉+𝑉𝜕𝑡^𝑚−𝑉_Na𝐺_Na, (4.3) with 𝐶_𝑚 denoting the membrane capacitance (F m⁻²). Since (4.2) has usually been solved under steady-state gating functions, which means that 𝑉≡𝑉(𝑥), then, passing to the limit as the thickness tends to zero, the steady-state conduction gives the characterisation 𝛽∼𝜎_𝑚membraneconductivity, (4.4) considering the Poisson equation to compute the electric field from the charge present into the system −∇⋅(𝜀∇𝑉)=𝐹𝑧𝑐. (4.5) 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, 𝜕𝑐𝜕𝑡+∇⋅𝐷∇𝑐+𝐹𝑧𝑘_𝐵𝑇𝑐∇𝑉=0, (4.6) and the adjoined Poisson equation (4.5), with 𝑘_𝐵 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 𝑐_𝑖 and 𝑐_𝑒 given by (3.26) and (3.24), respectively, 𝜕𝑐_𝑖𝜕𝑡=𝛼𝑐_𝑒−𝛽𝑐_𝑖, (4.7) then equality (3.28) reads ^∞_𝑗=0𝑑𝑝𝐷𝑓_𝑗+1𝑗+1+𝑞𝐷+𝛽(𝑡)𝑓_𝑗𝑥+ℎ+𝑝𝐷𝑡^𝑗=−𝛼(𝑡)𝛿_0,𝑖ℎexp−𝑏ℎ+𝜂−𝑎𝑥+𝜂²𝜎𝐷−^𝑖𝜀𝑡+()𝜉−𝑞𝐷exp−𝑏ℎ+𝜉−𝑎𝑥+²𝐷−𝜎_𝑒/𝜀−𝑞𝐷𝑡𝜉×𝜅₁𝛼(𝑡)1+𝐷2−ℎ𝜏exp(𝜏ℎ)+𝜅₂ℎ𝛼(𝑡)−𝐷2, (4.8) observing that 𝛼≡𝛼(𝑡) and 𝛽≡𝛽(𝑡). Therefore, it follows that 𝑀(𝑡):=𝑞/𝑝+𝛽(𝑡)/(𝑝𝐷),𝑃(𝑡)𝛼(𝑡)𝛿∶=−^0,𝑖ℎ𝑑𝑝𝐷exp𝑎−𝜂−𝑏ℎ+𝑎−𝜂𝑝+𝜂²−𝜎_𝑖,𝑁(𝑡)1𝜀𝐷−𝑞𝐷𝑡∶=)()𝑑𝑝𝐷𝜉exp(𝑎−𝜉−𝑏ℎ+𝑎−𝜉𝑝+𝜉²−𝜎_𝑒×𝜅𝜀𝐷−𝑞𝐷𝑡₁𝛼(𝑡)1+𝐷2−ℎ𝜏exp(𝜏ℎ)+𝜅₂ℎ𝛼(𝑡)−𝐷2. (4.9) 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 Section 4, 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.

Appendix

Following [42], we briefly present the main ideas of the theory for second-order linear differential equation with a regular singularity.

Consider the second-order linear differential equation 𝑦^′′+𝑝(𝑧)𝑦^′+𝑞(𝑧)𝑦=0. (A.1) We say that 𝑧=0 is a regular singular point for (A.1) if there exist 𝐹 and 𝐺 analytical functions at 𝑧=0 such that 𝑝(𝑧)=𝐹(𝑧)𝑧(𝑧)=𝐺(𝑧),𝑞𝑧². (A.2) We recall that 𝑢 is an analytic function at 𝑧=0 if 𝑢 is the sum of its Taylor series at 𝑧=0, in some neighbourhood of that point, that is𝑢(𝑧)=^∞_𝑛=0𝑢^(𝑛)(0)𝑧𝑛!^𝑛. (A.3) We now proceed with the following particular case (the Schrödinger equation for one particle, in spherical coordinates): 𝑑²𝜓𝑑𝑧²+1−𝜆−𝜆^′𝑧𝑑𝜓+𝑑𝑧−𝑘²+2𝛼𝑧+𝜆𝜆^′𝑧²𝜓=0, (A.4) where 𝜆, 𝜆^′, 𝛼, and 𝑘 are constants. Here, 𝐹(𝑧)=1−𝜆−𝜆^′,𝐺(𝑧)=−𝑘²𝑧²+2𝛼𝑧+𝜆𝜆^′. (A.5) For this equation 𝑧=0 is a regular singular point. Since 𝜆 and 𝜆′ are the solutions of the so-called indicial equation associated to (A.4), precisely 𝑟²+(𝐹(0)−1)𝑟+𝐺(0)=0, (A.6) then 𝜓₁=𝑧^𝜆𝑢₁(𝑧),𝜓₂=𝑧^𝜆′𝑢₂(𝑧) (A.7) are solutions for (A.4), where 𝑢₁ and 𝑢₂ are analytic functions at 𝑧=0. Setting 𝜓=𝑧^𝜆𝑓(𝑧), with 𝑓(𝑧)=𝑢₁(𝑧)+𝑧^{𝜆′−𝜆}𝑢₂(𝑧), then 𝑓 satisfies 𝑓^′′+1+𝜆−𝜆^′𝑧𝑓^′+2𝛼𝑧−𝑘²𝑓=0. (A.8) Now setting 𝑓(𝑧)=𝑒^−𝑘𝑧𝐹(𝑧), we obtain the following equation for 𝐹: 𝐹^′′+1+𝜆−𝜆^′−2𝑘𝑧𝐹^′−𝑘1+𝜆−𝜆^′−2𝛼𝑧𝐹=0, (A.9) which admits one analytical solution at 𝑧=0 (since its indicial equation is 𝑟(𝑟−𝜆′−𝜆)=0). Consequently, we get the solution 𝜓=𝑧^𝜆𝑒^−𝑘𝑧𝐹(𝑧) for (A.4). We obtain the explicit formula for the solution substituting the general series form 𝐹=∑^∞_𝑛=0𝐹_𝑛𝑧^𝑛 into (A.9). Setting 𝑤=𝑧/2𝑘, 𝑐=1+𝜆−𝜆′, and 𝑎=(1/2)(1+𝜆−𝜆′)−𝛼/𝑘, we can rewrite (A.9) into the following form: 𝑤𝑑²𝐹𝑑𝑤²+(𝑐−𝑤)𝑑𝐹𝑑𝑤−𝑎𝐹=0, (A.10) which is the so-called confluent hypergeometric equation.

Acknowledgments

Special thanks are due to Pedro C. Miranda, Carlos Farinha, and Ana Sebastião for some useful discussions on the subject and also to Catherine Pereira and Tiago Pereira who revised the english paper. The authors would like to thank the reviewers for their valuable suggestions. A. R. Domingos’s research was partially supported by Fundação para a Ciência e Tecnologia, Financiamento Base 2010, ISFL-1-209.