Improved Simulation of Electrodiffusion in the Node of Ranvier by Mesh Adaptation.

In neural structures with complex geometries, numerical resolution of the Poisson-Nernst-Planck (PNP) equations is necessary to accurately model electrodiffusion. This formalism allows one to describe ionic concentrations and the electric field (even away from the membrane) with arbitrary spatial and temporal resolution which is impossible to achieve with models relying on cable theory. However, solving the PNP equations on complex geometries involves handling intricate numerical difficulties related either to the spatial discretization, temporal discretization or the resolution of the linearized systems, often requiring large computational resources which have limited the use of this approach. In the present paper, we investigate the best ways to use the finite elements method (FEM) to solve the PNP equations on domains with discontinuous properties (such as occur at the membrane-cytoplasm interface). 1) Using a simple 2D geometry to allow comparison with analytical solution, we show that mesh adaptation is a very (if not the most) efficient way to obtain accurate solutions while limiting the computational efforts, 2) We use mesh adaptation in a 3D model of a node of Ranvier to reveal details of the solution which are nearly impossible to resolve with other modelling techniques. For instance, we exhibit a non linear distribution of the electric potential within the membrane due to the non uniform width of the myelin and investigate its impact on the spatial profile of the electric field in the Debye layer.

Introduction

Since the pioneer work of Hodgkin and Huxley [1], mathematical modelling of the electric activity of neurons has become an important tool to investigate the nervous system. While most models have relied on the cable theory formalism based on the strong analogy between neurons and electric circuits, limits of this formalism include its incapacity to account for fluctuations of ionic concentrations or to describe the electric field beyond the membrane. Extensions of the cable theory formalism have been developed [2, 3], notably in the NEURON software environment, in order to describe changes in ionic concentration as well as the impact of intracellular chemical reactions [4]. However, these extensions still fail to accurately describe the electric field beyond the membrane or the distribution of ionic concentrations in the Debye layer [5].

Solving the Poisson-Nernst-Planck (PNP) partial differential equations is a promising approach to overcome these limitations and model the evolution of ionic concentrations and of the electric field in neural structures such as axons, nodes of Ranvier, dendritic spines or the synaptic cleft ([5-9]). This strategy, not relying on oversimplifying assumptions such as the charge difference between intra and extracellular media being localized at the membrane, offers a potential spatial resolution in the nanometer range [7].

The FEM requires a spatial partition of the computational domain with a mesh of either tetrahedral or hexahedral elements and looks for a piecewise polynomial function on each of these elements which is the best approximation of the physical solution. In general, this renders the FEM unable to describe phenomena with spatial extents much smaller than the size of a mesh element which can potentially lead to exploding computational costs. For instance, the description of a volume of a cubic micron with an uniform resolution of ten nanometers would require a mesh of roughly one million elements. Modelling the evolution of four ionic concentrations as well as of the electric field on such a mesh using piecewise linear polynomials would imply having to solve 16 millions equations at each time step leading to either unreasonably slow or plainly infeasible calculations.

A clever way to reduce the computational cost without compromising on solution quality is to concentrate the mesh elements where the solution exhibits abrupt variations while tolerating a coarser mesh elsewhere. Such meshes, as the one used in [5], are created and tailored (adjusted) based on the knowledge of the behaviour of the solution, we will call them “tailored meshes”. It is also possible to rely on the mesh adaptation method which is based on error estimates derived from preliminary solutions [10-12]. Mesh adaptation is successfully applied in many industrial problems. While its first goal is to increase accuracy of the simulation for a given numerical effort, one can also view this approach as an automatic and objective method of constructing meshes based on intrinsic properties of the simulated quantities. As opposed to a “tailored mesh”, mesh adaptation doesn’t rely on the user a priori beliefs about the solution properties and remain easily applicable to complex geometries.

In neural structures, ionic concentrations as well as the electric potential are almost uniform over large parts of both the intracellular and extracellular space while experiencing steep variations in the vicinity of the membrane or in regions of high current density like the node of Ranvier or synapses. In the absence of transmembrane flux, variations in ionic concentrations and membrane potential are localized in the so called Debye layer [13] with a characteristic length of about 1 nm. These observations make the use of adapted meshes particularly promising in the field of computational neuroscience.

In the present paper, we first provide a detailed discussion on how to best apply the FEM to the resolution of the PNP equations in the field of computational neuroscience bearing in mind to make the use of our method by other research groups as straightforward as possible. Then, using a simple 2D geometry without transmembrane currents, we show that some kind of mesh adaptation (by an automatic process or by a “manual” positioning of the nodes) is necessary to obtain accurate solutions as resolution on an uniform mesh leads to sizeable numerical errors and for insufficiently refined meshes failed to converge to the exact value of the membrane potential obtained analytically. In the absence of transmembrane current, simulations based on automatic mesh adaptation presented here concentrate mesh nodes in the first nanometer from the membrane (in the Debye layer) exactly where the variations of the solution are expected to occur. We also compared the quality of solutions obtained on tailored meshes in which the node density decreased geometrically as a function of distance from the membrane (similar to the meshes in [5]). While both strategies led to important concentration of mesh nodes in the Debye layer, for a similar number of mesh nodes, the simulation on adapted mesh led to important gain in accuracy (in certain cases we observed a reduction by a factor of almost ten of the numerical error).

We finally model the generation of an action potential in a 3 dimensional model of a node of Ranvier. This establishes the feasibility of applying this approach to complex geometries thus to physiologically relevant investigations of structures such as dendritic spines, synaptic clefts or multiple nodes of Ranvier in nerve bundles [14]. Beyond that, our simulations reveal details of the solutions that cannot be captured by other approaches. For instance, we show that in the nodal region in the vicinity of the membrane, the electric potential can deviate by as much as 0.5 mV from neutrality (larger than predicted by models with uniform membrane width [5]) and that the electric potential varies in a highly non linear way in the membrane at the junction of the myelinated and non myelinated region. To our knowledge, our approach is the first to simultaneously describe the nanometric Debye layer and the complex geometry of a node of Ranvier.

It is our belief that those two tests provide a good illustration that simulations based on adapted meshes are the most efficient from a numerical point of view. Moreover, it can capture unsuspected minute details of the solution that cannot be detected without a priori knowledge. It is possible to argue that the use of adapted meshes could be neglected for oversimplified geometries as “manual” positioning of the nodes can also lead to good results. However, even in those simple cases mesh adaptation reduces the numerical efforts leading to shorter simulation times for a fixed accuracy level. More importantly, adapted meshes can be used on complex geometries for which “tailored meshes” are difficult or plainly impossible to construct.

Modelling methodology

Electro-diffusion theory and model

The electrodiffusion model solves the system of partial differential equations including the Nernst-Planck equation which describes the fluxes of n* ionic species subjected to electro-diffusion, as well as the Poisson equation which relates the electric potential V to the distribution of the electric charges. Here c, D and z are the concentration, diffusion coefficient and valence of the ionic specie k respectively. The coefficient α is defined by α = RT/Fz, where R is the perfect gas constant, F the Faraday constant, T the absolute temperature and ε is the dielectric constant of the medium Ω. If we define F, by as the flux of ionic specie k, then Eq (1) can be rewritten as the following continuity equation

We solved Eqs (1) and (2) on a computational domain Ω composed of three parts: the intracellular space, the membrane and the extracellular space. The diffusion coefficients of ionic species in the intra and extracellular spaces are taken to be the diffusion coefficients in water (see Table 1 and [15]) while these coefficients are set to zero on the membrane where the ionic concentrations themselves are also taken to be zero. The dielectric permittivity constant in both the intracellular and extracellular spaces is taken to be the dielectric permittivity constant of water ε0 ε where ε0 is the dielectric permittivity constant of vacuum and ε ≈ 80 is the relative dielectric permittivity of water [16] (see Table 1) while the dielectric permittivity of the membrane is given by ε0 εmem where εmem is the relative dielectric permittivity of the membrane which value will be discussed below.

Table 1

Electrodiffusion parameters.

R	8.31454 J ⋅ mole−1 ⋅ K−1	Perfect gas constant
F	96485 C ⋅ mole−1	Faraday constant
T	279.450 K	Absolute temperature
ε0	8.88541 ⋅ e−12 C ⋅ m−1 ⋅ V−1	Vacuum electric permittivity
εw	80	Water relative dielectric permittivity
εmem	40	Membrane relative dielectric permittivity
[K+]i	155mM	Initial intracellular K+ concentration
[Na+]i	12mM	Initial intracellular Na+ concentration
[A−]i	167.02mM	Initial intracellular anion concentration
[K+]o	4mM	Initial extracellular K+ concentration
[Na+]o	145mM	Initial extracellular Na+ concentration
[A−]o	149mM	Initial extracellular anion concentration
DK	1.96μm2 ⋅ ms−1	K+ diffusion coefficient
DNa	1.33μm2 ⋅ ms−1	Na+ diffusion coefficient
DA	2.00μm2 ⋅ ms−1	Anion diffusion coefficient.

While Eqs (1) and (2) could be used to describe an arbitrarily large number of ionic concentrations, we limited ourselves to a minimal set of cations K+, Na+ and a generic anion species A−. The first two are needed for the description of an action potential, whereas the last one is needed to obtain realistic net charge, osmolarity and cytosol resistivity [6]. Initial concentrations in the intracellular and extracellular domains are given in Table 1. For both concentrations and electric potential, we impose fixed Dirichlet conditions on the outermost boundary of the extracellular space and denote this boundary by . For ionic concentrations, we apply for each ionic species k = 1, ⋯, n*, the non-homogeneous boundary condition where take the initial value of the extracellular concentration of ionic specie k, whereas the voltage obeys to the following homogeneous Dirichlet condition On the outermost boundaries of the intracellular space denoted by and , we apply the following homogeneous Neumann boundary conditions for each ionic species k = 1, ⋯, n* and for the potential The boundary of the domain Ω is then given by . Moreover on the electrolyte-membrane interfaces, additional homogeneous Neumann conditions are applied for concentrations, except around the node of Ranvier where non zero flux is applied for concentrations. Denoting by the two cylindrical surfaces defining the intracellular-membrane and the extracellular-membrane interfaces at the node of Ranvier, we apply for every ionic specie where functions f are computed using the conductance of voltage-gated channels. The details of flux computation and the dynamics of the gating variables follow the standard Hodgkin-Huxley scheme and are given in the supplemental methods (S1 Methods). The vector n is the unit outer normal at the associated boundaries.

Weak formulation of the Electro-diffusion model

Multiplying Eq (4), for each ionic species k = 1, ⋯, n*, by proper test functions ψ and integrating over the domain Ω, we obtain For the electric potential, we multiply Eq (2) by the test function φ and obtain Applying integration by parts on Eqs (10) and (11) yields where . Taking into account Dirichlet boundary conditions Eqs (5) and (6), we denote by and functional spaces of concentrations and electric potential, respectively Reorganizing the previous system, we obtain the following so called weak formulation of the PNP equations: Find c ∈ H1(Ω) with , k = 1, ⋯, n*, and such that

Finite element approximations

The finite element treatment of the PNP equations first requires a spatial discretization of Eqs (13) and (14) by representing unknown functions c and V as well as test functions ψ and φ by piecewise polynomials on the mesh of simplicial ⊤ partitioning the domain Ω. In fact we consider the following finite approximation spaces and of piecewise polynomials of degree q and p defined relatively to the mesh : where is the space of polynomials of degree k* on ⊤. The discrete weak formulation in space of Eqs (13) and (14) is then given by: Find approximated solutions with and such that where is the interpolation of the initial value function onto .

For the discretization in time, we consider a second order backward difference formula (BDF2) (also known as Gear time stepping scheme, it is an implicit multistep marching scheme) [17]. This numerical procedure, for solving an ordinary differential equation of the type can be exemplified as follows: To obtain the approximation y of the solution y at time step t = t + Δt, solve the following equation which depends on values that are yet unknown. This implicit scheme is of second order accuracy, meaning that the error scales with . Applying it to the discrete weak formulation in space Eqs (17) and (18), we seek approximation of solutions c and V as follows: Given functions , , and , for times t and t respectively, find with and such that Expressing trial functions and as weighted sum of basis functions and for and respectively, and taking test functions as ψ = ψ, j = 1, …, N and φ = φ, j = 1, …, N, the system Eqs (21) and (22) is rewritten as an algebraic system of (n* × N)+N equations where and are vectors associated with the coefficients and in Eq (23). The coupling term in Eq (21) (third term on the left hand side) makes this an algebraic nonlinear system of generic form To solve such problem at time step r+1, the Newton-Raphson method is used. Starting with U1 = U (the solution at the preceding time step), a first order approximation of the linearized system is solved iteratively, and U is obtained at the end of this process (the loop is stopped when the correction δ is small enough). In our case, the main difficulty in the definition of the tangent matrix , arises from the coupling term. Introducing δ, δ and neglecting the second order terms, we get

Algorithm 1: Finite element resolution

1. Given a mesh, the initial time t0, , …, and V0 the initial conditions. Put r = 0;

2. While [final time is not reached]

• Compute a time step length Δt, t = t+Δt

• Compute the flux f by solving a system of ordinary differential equations (see S1 Methods)

• Put , V1 = V and s = 1

• While [desired tolerances on and δ are not reached]

- Construct the matrices , D

- Construct the vectors

- Solve the linear system

- Put

- Set s = s + 1

• End while

• Compute the approximation at time t

3. End while

Using this expression in Eqs (21) and (22) and the corrections δ, δ, the system Eq (24) is described by where, at time step r+1, we take as a starting state: and V1 = V, the matrices and vectors are defined respectively by α, β and γ are the coefficients of the finite difference Eq (20). The complete sequence of approximation is presented in Algorithm 1.

It is possible to have the length of the time step Δt changed as the calculation proceeds according to the dynamics of the system to speed up the simulation. Here we used the fact that the number of iterations for the resolution of one time step increases with Δt to reach a target range [Nmin, Nmax] of number of iterations for each time step. If Nit is the number of iterations at t, then for t we define

For the discretization of the system Eqs (13) and (14) we chose quadratic interpolation for all spaces, therefore k* = p = q = 2 and . Although the use of linear interpolations (i.e. ) could be viewed as less demanding (since it leads to smaller algebraic system), the higher degree of precision of quadratic interpolation makes it possible to use coarser meshes for the same precision.

The use of explicit time schemes for the system of Eqs (17) and (18) was excluded as it leads to conditional stability of the method and the Courant–Friedrichs–Lewy (CFL) condition [18] needed to insure stability would impose very small time steps. A simpler implicit time marching scheme could be used, for instance taking corresponds to the backward Euler scheme (BDF1). Once again, even if this could be viewed as a reduction of computation, the result would be opposite since the diminution in precision imposes the use of smaller time steps therefore an increase in the total number of time steps.

In the same manner, simplifications in the construction of the matrix involved in the left hand side of Eq (25) could be seen as reducing the assemblage, therefore decreasing computation time. However, this would lead to a degradation in the convergence rate of the iterative method thus to an increase in non linear iterations and overall calculation time. One such simplification could be to neglect the matrix contribution , still giving a converging method. Neglecting both terms of the linearisation ( and ) produces a diverging method since a large part of the information related to the coupling is ignored.

Mesh adaptation

Electric potential in neural structures is characterized by localized and fast transient behaviour. In such cases the accuracy of the numerical approximations often deteriorates due to phenomena such as: local singularities (like those arising from re-entrant corners of domains), excessively small zones of heterogeneous material, presence of boundary layers or sharp moving fronts. An obvious strategy to improve the quality of the solution is to refine the resolution grid but in situations where different spatial scales are involved, refining in a global manner will lead to an excessive work load. Mesh adaptation improves the quality of FEM solutions thanks to local operations such as edge refinement (Fig 1A), node elimination (Fig 1B), edge swapping (Fig 1C) and node displacement (Fig 1D).

Fig 1

a) Edge refinement. b) Node elimination. c) Edge swapping in 2D. d) Node displacement.

We will briefly detail the main ideas behind our hierarchical error estimator and our adaptive remeshing strategy and we refer to [11, 12] for more details. Let denote a triangulation of the domain Ω and ⊤ its elements. Let u be the exact solution, usually unknown and suppose that we have computed a finite element approximation of degree q of the solution denoted . We now make the assumption that we can build, starting from , a new approximation of degree q + 1 that is at least slightly more accurate in the sense that with β < 1. The construction of is therefore crucial and will be described later on. For the moment, we assume its existence and the triangle inequality gives thus and the error can then be controlled by the right-hand side in some appropriate norm (L2 or H1-norm for instance).

Now the idea is to define such as: where can be seen as a degree q + 1 correction to the degree q approximate solution , thus the name hierarchical error estimator. The error bound Eq (29) reduces to

In order to construct our method requires an accurate approximation of the gradient (∇u) of the solution (denoted ), which we also require to belong to the space of piecewise continuous functions whose restriction to any element ⊤ of belongs to the space ). For this purpose, there exists gradient recovery methods producing accurate approximations of ∇u. We will use the one described in [19].

Then can be easily computed using and the recovered gradient . Indeed, to enrich in order to get an approximation of degree q + 1, we need, as in a Taylor expansion, its derivatives of order q + 1. Obviously, the (q + 1)th derivatives of vanish but we can use the qth derivatives of the recovered gradient . This means that the (q + 1)th derivative of should coincide with the appropriate qth derivative of .

In the general case, a well posed linear system can be built on each element, whose solution completely defines c(. Now, if we can control the correction c(, then from Eq (31) we can also control the error and this will serve for the construction of an adapted mesh.

The first objective is to reach a prescribed global target level of error (eΩ) in L2-norm. Since we are modifying the mesh locally, this objective will be reached also locally by imposing on each element This way, the error is distributed proportionally to the element area (or volume in 3D) and if this target value can be reached on each element, then the global error eΩ will be attained. These local target errors are imposed using edge refinement and node elimination (Fig 1A and 1B). For edge refinement for example, we build for each edge e of the mesh the patch of elements containing e (Fig 1A). On that patch, the local target error should be and the estimated error is We then cut that edge by adding a new mid-side node and we reinterpolate the recovered gradient at the new node which is necessary for the computation of c( on the newly created elements. We then compute the error estimate on the new patch and choose, between the initial patch and the new one, the one for which the estimated error is closest to the target Eq (32). A similar procedure is used for node elimination and more details are given in [11, 12].

We also want to reach the global target error value with a minimum number of nodes. As a second objective, we therefore try to achieve some form of equidistribution of the error by minimizing, as in [20], the H1-seminorm of the error which is approximated by the quantity: This minimization is also performed locally but now using edge swapping and node displacement. To move node for example, we construct the patch of elements sharing that node and determine its new position in order to minimize the quantity A gradient method can be used for the minimization. The node is only allowed to move inside the patch to prevent element inversion. A similar procedure is used for edge swapping. Remarkably, the minimum on the H1-seminorm of the error cannot be achieved without somehow reorienting and stretching the elements in appropriate directions, therefore leading to anisotropic meshes as shown in [12, 20].

In summary, mesh modification is done by sweeping the nodes (for node elimination and node displacement) and the edges (for edge division and edge swapping) a number of times until the two above objectives are approximately satisfied.

Algorithm 2 summarizes the method used for mesh adaptation.

Algorithm 2: Mesh adaptation method

1: Given an initial mesh .

2: While [mesh variations occurs]

3:  Solve to obtain a solution u on mesh .

4:  Compute an estimation of the error Eq (31).

5:  Adapt the mesh using local operations to obtain .

• Reach a prescribed global target level of error in L2-norm thanks to edge refinement and node elimination.

• Minimize the H1-seminorm of the error using node displacement and edge swapping (error equidistribution).

6:  Set i = i+1.

7: End while

Algorithm 2 was performed on an initial mesh generated from COMSOL Multiphysics, a commercial software.

Finally, the method was implemented using the finite element library MEF++ developed at the Groupe Interdisciplinaire de Recherche en Éléments Finis (GIREF) (see http://giref.ulaval.ca/mef.html).

Results

Gain in accuracy through solution based mesh adaptation

While the FEM and the finite volume method have been used to describe various problems in neuroscience ([5, 6]), to the best of our knowledge no attempt has been made to quantify the relationship between the numerical error and the mesh size in this context. For the sake of simplicity and in order to be able to compare our results with exact solutions, we first considered the simplest of situations, a two-dimensional geometry with linear membranes and no transmembrane flux (Fig 2A). In this scenario, variations in ionic concentrations and electric potential are localized near the membrane where the electric potential is well described by the Poisson Boltzmann equation with k ≈ 1.38 × 10−23 J ⋅ K−1 the Boltzmann constant and N ≈ 6.022mol−1 the Avogadro number. In this context denotes the concentration of ionic specie k at infinity (far from the membrane) taken to be equal to its initial value. Linearizing Eq (33) and assuming electroneutrality away from the membrane, one obtains Solving Eq (34) yields the characteristic length of the region in which the electric potential and ionic concentrations vary steeply, i.e. the Debye layer (λ), according to [13] For the specific values of our problem this gives a Debye length of approximately 1.1 nm.

Fig 2

a) Schematic of the two dimensional model used to compute the numerical error of the method. The membrane potential was computed by the difference of electric potential at each side membrane (see the two circles). b) Example of a computational grid obtained with the mesh adaptation method. c) Distribution of mesh nodes in the intracellular space as a function of the distance from the membrane. Results obtained for the iterations Algorithm 2 as well as for a tailored mesh. d) Membrane potential error (see Eq (37)) on different meshes. The error was computed as a function of the number of computation nodes in the mesh for different strategies: uniform meshes with a uniform refinement, tailored meshes with increased node density near the interface of the membrane and intra (extra) cellular space and meshes obtained through mesh adaptation method.

These observations about the Debye layer led us to formulate three hypotheses: 1) The mesh adaptation method should concentrate a large proportion the mesh elements within 1 nm of the interfaces between the membrane and the intra-extra cellular spaces, 2) Mesh node density should decrease in a way that is approximately exponential with distance from the membrane, 3) Given that variations occur only in the direction perpendicular to the membrane, the adapted mesh should be very anisotropic.

To test these hypotheses, we iteratively applied the mesh adaptation method and solved the PNP equations on the domain (see Algorithm 2). This process led to a very anisotropic mesh with nodes concentrated near the interfaces between membrane and intra-extra cellular spaces (Fig 2B). We created a tailored mesh (see details below) as an alternative to the mesh adaptation method. In Fig 2C we present the nodal distribution, in the intracellular space, as a function of distance from the membrane for the meshes produced by the adaptation process and our tailored mesh. This figure also depicts the evolution of the adaptation process. Here the process was stabilized (i.e. no notable changes in mesh topology observed) after 15 iterations. This clearly demonstrates (Fig 2C) the densification near the membrane in a zone less then 0.25 nm from the membrane. Notice the overshoot of the sixth mesh with a peak of nearly 55% of the nodes, which is corrected in the next iterations by redistributing the superfluous nodes in the nearby zone, leading to a mesh with around 20% of the nodes in the first 0.25 nm and almost 45% of the nodes in the first nm.

These results can be interpreted as the mesh adaptation algorithm detecting the necessity for the size of some mesh elements to be smaller than 1 nm. However, for a regular mesh to have elements of this size in a 2D geometry of 1 μm2 surface area would require 106 elements bearing a prohibitive computational cost. A tailored mesh as proposed by J. Pods et al. [5] is more efficient compared to a uniform mesh. However it requires an a priori knowledge of the solution and in the case of complex geometries, such an approach can become impracticable since properties of the solution may be difficult, if not impossible, to infer in advance. As for the adaptive method, in any cases, it will automatically densify meshes on sensitive areas with respect to the accuracy of the numerical solution. As an indication we added, in (Fig 2C), the nodal distribution for a tailored mesh which is similar to the mesh proposed in [5]. We created this mesh by manually specifying regions of high node density near the membrane and regions of low density elsewhere. The mesh was constructed so that the node density decreases in a geometric manner as a function of the distance from the membrane.

To compare the numerical errors of solutions obtained with different meshing strategies, we took advantage of the easily computed value of the membrane potential where L is the length of the membranes, dintra the distance between the two membranes, dmem is the width of the membrane (see Fig 2A). The electric permittivity of the membrane is ε and Q is the mean electric charge density in the intracellular medium given by We compared this exact value to the membrane potential computed by the FEM approach at the center of the membrane. More precisely we computed the absolute difference between the variation of Vmem and the variation of our approximation V on each side of the membrane at control points and (see Fig 2A):

The rate (order) of convergence, with respect to the time step and the spatial size of the elements are relatively standard theoretical results (see [21] for instance). The FEM (for continuous coefficients and a regular geometry), would give an error of order k⋆ with respect to space since we are using degree k⋆ polynomial approximation in space and of order 2 in time since we are using an order 2 marching scheme (here BDF2). Because we are not fulfilling all the conditions here (our coefficients are discontinuous), and since we are using quadratic polynomial approximation, as presented by Ying and Benzhuo [22] (see also [23]), we “lose” an order of convergence in space leading to a rate of convergence of 1 in space. It should be noted that using a linear approximation in space would lead to a rate of converge in space of 1/2, which strengthens the argument in favour of the use of quadratic interpolation. The approach of mesh adaptation has no effect on the convergence rate since the theoretical results don’t require the discretization parameters to be constant. Of course, in the case of non uniform meshes, error graphs cannot be based on the size of elements which varies throughout the domain, but rather on the number of nodes, which becomes the comparing quantity. Since we have a rate of convergence of 1, the error is expected to be inversely proportional to the number of nodes. We also illustrated the spatial rate of convergence of steady-state solutions (Fig 2D) as we have slopes of 1 for the adapted mesh method as well as for the tailored meshes, indicating that both approaches can be used with confidence in this setting. Notice that the curve for the uniform meshes cannot give us this information as excessively refined meshes would be needed to exhibit a slope of 1. Indeed, the values of membrane potential computed on regular meshes obtained by successive global refinements failed to converge to the exact value one for a number of elements up to 105.

Contrastingly, when we applied the mesh adaptation technique to meshes of different sizes, the error was less 10−4 mV for as little as 695 mesh elements (Fig 2D). Using different tailored meshes (obtained by increasing node density in each region) the approximation of the membrane potential converged to the theoretical value with the error reaching 10−4 mV for 1458 elements (Fig 2D) therefore needing more then two times the number of nodes used by the adapted mesh for the same order of accuracy.

Since we have a convergence rate of 1 in both cases, the gap in accuracy and in number of nodes will remain, even if we refine the meshes (which corresponds to increasing the number of nodes). Using our initial tailored mesh as a reference, the mesh adaptation method reduces the number of nodes for a given precision by a factor nearly equal to two (Fig 2D). For an accuracy of about 10−5, the adaptive method requires approximately 5000 nodes while we would need to refine our initial tailored mesh until it has more than 10000 nodes.

Using tailored meshes with better positioning and density of nodes could produce accuracy arbitrarily close the one obtained using adapted meshes, which are in a sense optimal as discussed above. However, this raises the question of how to determine, with the same efficiency and precision, the zones and the densities of nodes in those regions. Pretending to establish optimal targeted geometrical zones, densities and to control the total number of nodes without tools comparable to those used by the adaptation method (i.e. without a detailed knowledge of the behavior of the error) seems inconceivable.

Application of mesh adaptation in a model of action potential generation in a node of Ranvier

The advantages of the mesh adaptation method being demonstrated in the simplistic instance of a 2D geometry, we now apply this strategy to the simulation of an action potential in a node of Ranvier. The details of the geometry and most parameters of the model were taken as in [6] while the leak conductance was taken as in [5], see also Fig 3A and Table 2. The action potential was provoked by imposing a [Na+] flux in the nodal region during the first 0.5 ms of the simulation. The time course of membrane potential and of intracellular ionic concentrations are shown in Fig 3B–3D and are comparable to results obtained by Lopreore et al. [6].

Fig 3

a) Schematic of the three dimensional geometry used in the model of a node of Ranvier. b-d) Time of electric potential (b), potassium concentration (c) and sodium concentration (c) during an action potential at the spatial point indicated by a circle in a.

Table 2

Additional parameters for the node of Ranvier.

d	0.434μm	Radius of the axon
Dn	0.02μm	Thickness of the unmyelinated part of the membrane
Dm	0.406μm	Thickness of the myelinated part of the membrane
L	4μm	Length of the axon section
Lnode	0.7μm	Length of node in which currents are applied
gKL	0.435mS ⋅ cm−2	Conductance density of leak K+ channels
gNaL	0.065mS ⋅ cm−2	Conductance density of leak Na+ channels
g¯Kv	36mS ⋅ cm−2	Conductance density of voltage-gated K+ channels
g¯Nav	120mS ⋅ cm−2	Conductance density of voltage-gate Na+ channels
dur	0.5ms	Duration of the stimulus
Dur	10ms	Duration of the simulation

We adjusted the relative dielectric permittivity of the membrane to obtain an equivalent electric capacitance of 2μF/cm2 in the node and 0.5μF/cm2 in the myelinated part of the membrane as in [6]. We made the conversion according to the formula giving electrical capacitance of a cylindrical shell where d is the radius of the axon and D the width of the membrane. This gave the unphysiological value of εmem = 40 (as opposed to the experimentally measured value between 2 and 10 [24]). This might be due to the fact that the width of membrane in the nodal region is larger than the physiological value (20 nm vs ∼ 5–10 nm [25]). The specific values of these parameters have however no impact on the general results of the present paper regarding the benefits of the mesh adaptation method.

Highlighting regions of interest of the solution, the mesh adaptation method concentrates mesh elements in and near the unmyelinated part of the membrane while the elements outside of this region are very anisotropic and especially scarce in the extracellular space (Fig 4A). Beyond increasing the accuracy of the results (as described in the previous section), by concentrating the nodes of the grid in regions of interest, adapted meshes can reveal features of the solutions that would otherwise be missed. In particular, the frequent assumption that the electric field is constant in a cross section of the membrane (as is used to derive the GHK flux equation) is no longer valid when the width of the membrane is not uniform as in the present case (Fig 4B). Remark that solving the Poisson equation on complex geometries, especially in the presence of reentering corners, is a problem related to intrinsic mathematical difficulties so it would be difficult to obtain a good approximation of the solution by simpler methods.

Fig 4

a) Illustration of the three dimensional mesh obtained for the description of node of Ranvier by the mesh adaptation method. b) Electric potential on the whole domain taken at the peak (most depolarized time point) of the action potential illustrating its non linear distribution on cross-sections of the membrane. c-e) Enlarged view of the electric potential in small regions near the membrane illustrating that the model is able to describe the so-called Debye layer.

Eq (34) suggests that the amplitude of the deviation from neutrality of the electric potential in the Debye layer of the extracellular space is directly proportional to its gradient on the membrane (though the actual picture can be complicated by the presence of transmembrane ionic fluxes). This predicts a larger deviation near the corner between the myelinated and unmyelinated section of the membrane (Fig 4D) than near the middle of the myelinated section (Fig 4E).

Discussion

Solving the PNP equations with the FEM is a promising approach bound to gain wider use to describe electro-diffusion in neural structures [7]. This method is particularly useful to describe structures with complex geometries such as nodes of Ranvier, dendritic spines or the synaptic cleft and its surrounding [5, 6, 8, 9]. In such contexts, the FEM describes phenomena that other computational approaches fail to capture such as the fine scale distribution of ion concentrations and of the electric potential or the way their variations extend beyond the membrane.

However, discontinuities in dielectric permittivity and ionic diffusion coefficients occurring at the interfaces between the membrane and the intra-extra cellular spaces pose numerical challenges, as the solutions typically exhibit strong variations localized at the vicinity of these interfaces. As demonstrated here, this feature of problems naturally occurring in the field of neuroscience can lead to significant numerical errors if not addressed properly. The efficiency of “tailored meshes” have been demonstrated in the context of an axisymmetric model of propagation of an action potential in an axon with membrane of uniform width [5]. In the present paper, we show that the technique of the mesh adaptation method, which has a proven track record in the field of industrial mathematics, can be applied to describe electrodiffusion by solving the PNP equations on complex three dimensional geometries in neuroscience. In addition of yielding higher performance, this approach has the advantage of not requiring the user to manually specify mesh properties based on possibly incomplete a priori beliefs on the solution.

A potential limitation of the mesh adaptation method is that since this procedure is performed using a solution at a given time, the adapted mesh may no longer be optimal as the features of the solution change over time. In the instance of the propagation of an action potential along an axon for example, an ideal strategy would make the concentration of mesh elements follow the propagation of this action potential since this is where the solution varies most abruptly and where it is most important to describe the details of the solution. A possible way to achieve such a time evolving mesh would be to readapt the calculation grid at several time points of the simulation (at each tenth of millisecond for instance) as we hope to achieve in future research. Another potential limitation of using the FEM to solve the PNP equations is that the use of partial derivative equations implicitly assumes that ionic concentrations are continuous quantities and thus fails to describe the granularity arising from the discrete distribution of ions and of transmembrane proteins such as voltage-gated channels. It has however been demonstrated in [6] that stochastic equations can be coupled with partial derivative equations and capture the variance in the solution due to the stochastic behaviour of channels. It might of interest to go one step further and use a system of stochastic partial derivative equations to account for the probabilistic movement of ions. This could become important for ionic species with low concentration such as Ca2+ where the small number of ions in each spatial elements can violate the continuity condition on which the use of diffusion equations is founded.

Finally, the FEM applied to structures with geometry evolving over time leading to mesh deformation is relatively common (for instance, in biology, problems such as the dynamics of visicle shape transformations or mitral valve simulation, in engineering, the extreme deformation and buckling of structures). However the integration of the adaptation techniques for such problems is more recent (see [26]). Such techniques could be applied in the field of computational neurosciences, to describe the birth and growth of dendritic spines [27] or the swelling of cells [28], potentially helping to better understand these important phenomena.

Conclusion

In the field of neuroscience, experimental data is discovered at an ever increasing pace. Integrating this knowledge in comprehensive models is a great challenge promising to help us better understand the function of neural structures. In order for this to be achieved, the steep computational requirements of modelling micrometric structures with a spatial resolution sufficient to fully explain their function (often required to be in the nanometric range) will have to be addressed. While the raw power of parallel calculations performed on super computers will be required, we believe that many numerical methods developed in the field of industrial mathematics could allow to push back the numerical wall arising in computational neuroscience. As demonstrated in the present paper, the use of the mesh adaptation method greatly reduces the numerical cost of solving the PNP equations to model electro-diffusion while increasing the precision of the solution. We hope that this research will open the door for the modelling of larger structures such as axons spanning of several nodes of Ranvier or multiple axons in a nerve bundle. This could provide insights on how loss of myelin can affect signal propagation [14, 29].

Detailed information on the dynamics of voltage-gated channels and on the simulation of the generation of an action potential.

(PDF)

Click here for additional data file.