#### 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^{−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. 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
for some (see Figure 1). Let be a neighbourhood of the membrane defined as
for some . Thus and denote the intracellular space and the membrane surface, respectively. Define the external boundary
Let . The problem under study is defined by the system of parabolic equations (for details see [23])
where is the sodium diffusion coefficient ( m^{2} s^{−1}), [27]) and is the sodium permittivity ( F m^{−1}). 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^{−3}. The electrical conductivity , with , is considered homogeneous and constant at each subdomain (extra- and intracellular domains), measured in S m^{−1} [28]:
Other values for the diffusion coefficient and the electrical conductivity can be found in [29].

**(a)**

**(b)**

We assume the following boundary conditions: 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 , 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 the following condition holds: where . 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 with , for every (see Figure 2). The angular position 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).

The initial condition is assumed constituted by an averaged condition:
with the values of and depending on the physiological data (e.g., in a cat motoneuron mol m^{−3} and 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 and on
where represents the depolarisation voltage, is the Faraday constant ( C mol^{−1}, [3]), denotes the universal gas constant ( J (mol·K)^{−1}), 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: with . 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, , , , , and (see Section 4). Thus, we get a solution (see Section 3.1, for details): with given by (3.25), and where , , , 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 Consequently, from substituting this in (3.1) we obtain where and , for . Then there exists such that

*Remark 3.1. *The second-order ordinary differential equation for has the following solutions.

If , then there exist real constants , such that
if , then there exist real constants , such that ; if , then there exist real constants , such that
We observe that the validation of the phenomenal data yields the choice of 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 with nonnegative real constant. We will now analyse the second equation from (3.6), that is, Inspired by the theory for the confluent hypergeometric equations (see the appendix), we look for a function of the form that satisfies (3.10), where and , , , and are real constants to be chosen later. Solving (3.10) for , we get Thus, we choose (see [34]) with , and in order to get positive solutions. Solving (3.10) for , we obtain

###### 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 Applying the Taylor formula in , with , and denoting , we have Since there is no dependence in time, we get

Denoting and , it follows that From (3.15) it follows that Thus, introducing (3.19) into (3.20) and using (3.18), after some calculations we obtain , with Thus, and consequently

Notice that it still remains to find , , , , , and . Using (2.8), (3.9), and (3.23), we obtain (hence, ), , , and , concluding that

Finally, the initial condition (2.10) yields 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): with and . Here, we consider new unknown constants , , , , , , , 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, Next, we verify (2.7). Introducing (3.24) and (3.26) in (2.7), we obtain 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, where and

Introducing (3.29) into (3.27), we conclude that , where 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 These mean that the constants and are well defined, and and for all .

Notice that the constants , and are well defined observing that the constants , , 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 This yields the choice for . Finally, from (2.11) it follows that for all . Hence, the algebraic system for and is 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 are determined in (3.37) and (3.38), which come from the Boltzmann principle (2.11),(v) and 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 , 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 , a turning-on process , and an inactivating process , where and denote the activation and the inactivation of the sodium current, respectively [3, 4, 7]:
Their dynamics are described by the ODE system:
where , with and 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
with denoting the membrane capacitance (F m^{−2}). 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
considering the Poisson equation to compute the electric field from the charge present into the system
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, 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, then equality (3.28) reads observing that and . Therefore, it follows that , 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
We say that is a * regular singular point* for (A.1) if there exist and analytical functions at such that
We recall that is an analytic function at if is the sum of its Taylor series at , in some neighbourhood of that point, that is
We now proceed with the following particular case (the Schrödinger equation for one particle, in spherical coordinates):
where , , , and are constants. Here,
For this equation is a regular singular point. Since and are the solutions of the so-called * indicial equation* associated to (A.4), precisely
then
are solutions for (A.4), where and are analytic functions at . Setting , with , then satisfies
Now setting , we obtain the following equation for :
which admits one analytical solution at (since its indicial equation is ). Consequently, we get the solution for (A.4). We obtain the explicit formula for the solution substituting the general series form into (A.9). Setting , , and , we can rewrite (A.9) into the following form:
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.