Sensitivity of the electrical response of a node of Ranvier model to alterations of the myelin sheath geometry

Abstract

Nodes of Ranvier play critical roles in the generation and transmission of action potentials. Alterations in node properties during pathology and/or development are known to affect the speed and quality of electrical transmission. From a modelling standpoint, nodes of Ranvier are often described by systems of ordinary differential equations neglecting or greatly simplifying their geometric structure. These approaches fail to accurately describe how fine scale alteration in the node geometry or in myelin thickness in the paranode region will impact action potential generation and transmission. Here, we rely on a finite element approximation to describe the three dimensional geometry of a node of Ranvier. With this, we are able to investigate how sensitive is the electrical response to alterations in the myelin sheath and paranode geometry. We could in particular investigate irregular loss of myelin, which might be more physiologically relevant than the uniform loss often described through simpler modelling approaches.

Appendices

A Discrete version of the Poisson Nernst–Planck equations

We denote by $L^{2}(\Omega )$ the space of square integrable functions defined on $\Omega$. The Sobolev spaces $H^{m}(\Omega )$, $m \ge 0$, are the spaces of functions in $L^{2}(\Omega )$ with generalized partial derivatives belonging to $L^{2}(\Omega )$ up to order m (see Adams and Fournier 2003). The functional space $U_V$ containing the ionic concentrations and electric potential is defined as

$$\begin{aligned} V = \left\{ v \in H^{1}(\Omega ): v = 0 \text{ on } \Gamma _{D}\right\} ,\quad U_V = L^{2}(0, T_F; V)\cap L^{\infty }(0, T_F; L^{\infty }(\Omega )). \end{aligned}$$

1.1 A.1 Weak formulation of the Poisson Nernst–Planck equations

Multiplying the system of Eq. (1) by test functions $v, w \in H^{1}(\Omega )$ and integrating over the domain $\Omega$, we obtain the following weak formulation

$$\begin{aligned} \left\{ \begin{aligned}&\int _{\Omega }\left( \frac{\partial c_{i}}{\partial t} + \nabla \cdot {\textbf{J}}_{i}(c_{i},\phi )\right) vdx = 0,\;\;i=1,\ldots ,n_c, \\&- \int _{\Omega }\left( \nabla \cdot (\varepsilon \nabla \phi ) + F \sum _{i = 1}^{n_{c}} z_{i}c_{i}\right) wdx = 0. \end{aligned} \right. \end{aligned}$$

Integrating by parts and taking into account (2), (4) and (5) we have $\forall v,w \in V$,

$$\begin{aligned} \left\{ \begin{aligned}&\int _{\Omega }\frac{\partial c_{i}}{\partial t}vdx -\int _{\Omega }{\textbf{J}}_i(c_{i},\phi ) \cdot \nabla vdx \\&\quad = - \sum ^{n_{sd}}_{k=1}\int _{\partial \Omega _{k}}\left( {\textbf{J}}_{i}(c_{i},\phi )\cdot {\textbf{n}}_{k}\right) vds,\;\; i=1,\ldots ,n_c, \\&\int _{\Omega }\varepsilon \nabla \phi \cdot \nabla wdx -F\sum _{i = 1}^{n_{c}}\int _{\Omega }z_{i}c_{i}wdx = \sum ^{n_{sd}}_{k=1}\int _{\partial \Omega _{k}}\left( \varepsilon \nabla \phi \cdot {\textbf{n}}_{k}\right) wds. \end{aligned} \right. \end{aligned}$$

Since

$$\begin{aligned} \partial \Omega _k = \left( \partial \Omega \cap \partial \Omega _{k}\right) \cup \left( \partial \Omega _k\cap \Gamma \right) \quad \text{ and }\quad \left( \partial \Omega \cap \partial \Omega _{k}\right) \cap \left( \partial \Omega _k\cap \Gamma \right) = \emptyset \end{aligned}$$

taking into account (6) and introducing ${\tilde{D}}_{i} = D_{i}/\lambda _{i}$ we obtain the weak problem: find $c_{i} \in U_V$, $i = 1,\ldots , n_c$ and $\phi \in U_V$ such that

$$\begin{aligned} \left\{ \begin{aligned}&\int _{\Omega }\frac{\partial c_{i}}{\partial t}vdx + \int _{\Omega }(D_{i}\nabla c_{i} + {\tilde{D}}_{i}c_{i}\nabla \phi )\cdot \nabla vdx \\ {}&\quad = - \int _{\Gamma }f_{i}vds,\; \forall \, v \in V,\;\; i=1,\ldots , n_c,\\&\int _{\Omega }\varepsilon \nabla \phi \cdot \nabla wdx -F\sum _{i = 1}^{n_{c}}\int _{\Omega }z_{i}c_{i}wdx = 0,\;\forall \, w \in V. \end{aligned} \right. \end{aligned}$$

(9)

1.2 A.2 Numerical approximations

In order to numerically approximate solutions of (9) we introduce a time discrete formulation based on a finite difference time scheme and combine it with the finite element method for the spatial approximation. We consider a second order backward difference formula with, for convenience and simplicity, a constant time-step $\Delta t$. Denoting $t^n = n\Delta t$ and $\psi ^n = \psi (\cdot ,t^n)$, the time derivative of $\psi$ at time $t^{n+1}$ is approximated by

$$\begin{aligned} \frac{\partial \psi (\cdot ,t^{n+1})}{\partial t} \approx \frac{3\psi ^{n+1}-4\psi ^n+\psi ^{n-1}}{2\Delta t},\qquad n\ge 1. \end{aligned}$$

The first timestep being calculated using a backward Euler approximation.

For the spatial approximation, the finite element method is used. First, we introduce a mesh ${\mathcal {T}}_{h} = \{\top \}$ of simplicial $\top$ partitioning $\Omega$:

$$\begin{aligned} {\overline{\Omega }}_h = \bigcup _{\top \in {\mathcal {T}}_{h}}\top \,\subseteq {\overline{\Omega }}, \end{aligned}$$

with the usual restrictions on ${\mathcal {T}}_{h}$ (see Bathe 1996 or Ciarlet and Luneville 2009 for details). From the family ${\mathcal {T}}_{h}$ of triangulations of the domain $\Omega$ (indexed by h), we construct the family of finite dimensional vector spaces $V_{h}$. This set of standard continuous and piecewise functions is defined by

$$\begin{aligned} V_{h} = \left\{ v_h \in {\mathscr {C}}^{0}\left( \overline{\Omega _h}\right) ;\, {v_h}_{|\top } \in {\mathscr {P}}_{k}\left( \top \right) , \, \forall \, \top \in {\mathcal {T}}_{h}, \, v_h = 0 \, \text{ on } \Gamma _{D} \right\} \end{aligned}$$

where ${\mathscr {P}}_{k}(\top )$ is the space of polynomials of degree less or equal to k on $\top$. An index "h" is used to denote a spatial approximation in $V_h$. An approximation, in $V_h$, of $\psi$ at time $t^n$ will be noted $\psi ^n_h$.

Since an implicit time scheme is used for an approximation at time $t^{n+1}$, all terms (apart form the time derivatives) will be expressed at time $t^{n+1}$. Applied on (9), this results in sequences of strongly coupled (and possibly non linear) problems characterizing time-discrete unknown (ionic concentrations, electric potential, etc.). Details for the resolution of these non linear systems and the resulting algorithm, are presented in Dione et al. (2016).

1.3 A.3 Discrete formulation of the PN–P equations

To simplify the notation, let

$$\begin{aligned} f_{ih}^{n+1} = f_{i}(t^{n+1}, c^{n+1}_{ih}, \phi ^{n+1}_{h}). \end{aligned}$$

The implicit time scheme and finite element method applied on (9) produce a sequence of problems: for $n \ge 1$, given $c^{n-1}_{ih}$, $c^{n}_{ih}$,$\;\;i=1,\ldots ,n_c$, $\phi ^{n-1}_{h}$ and $\phi ^{n}_{h}$, find $c^{n+1}_{ih} \in V_{h}$ and $\phi ^{n+1}_{h} \in V_{h}$ such that

$$\begin{aligned} \left\{ \begin{aligned}&\int _{\Omega _{h}} \frac{3c^{n+1}_{ih} - 4c^{n}_{ih} + c^{n-1}_{ih}}{2\Delta t_n}v_{h}dx + \int _{\Omega _{h}}(D_{i}\nabla c^{n+1}_{ih}+{\tilde{D}}_{i}c^{n+1}_{ih}\nabla \phi ^{n+1}_{h})\cdot \nabla v_{h}dx \quad \\&\quad = - \int _{\Gamma _h}f_{ih}^{n+1}v_{h}ds, \; \forall \ v_{h} \in V_{h},\;\;i=1,\ldots ,n_c, \\&\int _{\Omega _{h}}\varepsilon \nabla \phi ^{n+1}_{h}\cdot \nabla w_{h}dx -F\sum _{k = 1}^{n_{c}} \int _{\Omega _{h}}z_{i}c^{n+1}_{ih}w_{h}dx = 0,\forall \; w_{h} \in V_{h}, \end{aligned} \right. \end{aligned}$$

(10)

where $\Gamma _h$ is the discrete version of the internal interface $\Gamma$. These systems are non linear, since they contain ${\tilde{D}}_{i}c^{n+1}_{ih}\nabla \phi ^{n+1}_{h}$, moreover, the right-hand side $f_{ih}^{n+1}$, a priori, depends on $c_{ih}^{n+1}$. We use a Newton–Raphson method for solving (10) (which amount to a first order Taylor development of the non linear terms), the reader is refered to Dione et al. (2016) for the complete details.

B Voltage-gated channels

The functions $n(t, \phi )$, $m(t, \phi )$ and $h(t, \phi )$ are time and voltage-dependent gating functions, taking values in [0, 1] defined by the following ODEs

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{dn}{dt} = \alpha _n(\phi )(1-n) - \beta _{n}(\phi )n, &{} \\ \displaystyle \frac{dm}{dt} = \alpha _m(\phi )(1-m) - \beta _m(\phi )m, &{} \\ \displaystyle \frac{dh}{dt} = \alpha _h(\phi )(1-h) - \beta _h(\phi )h, &{} \end{array}\right. } \end{aligned}$$

(11)

where, denoting $\phi _r$ the resting potential, the rate functions $\alpha _{(\cdot )}$ and $\beta _{(\cdot )}$ are defined by

$$\begin{aligned}{} & {} \phi _{r}=-55\; mV,\quad \gamma _v(\phi ) = \left( v + \phi _r - [\phi ]\right) /10 \\{} & {} \beta _n(\phi ) = \displaystyle 0.125 e^{\gamma _0(\phi )/8}, \quad \alpha _n(\phi ) = {\left\{ \begin{array}{ll} 0.1\displaystyle \gamma _{10}(\phi )\left( e^{\gamma _{10}(\phi )} - 1\right) ^{-1} &{} \text{ if } \; \displaystyle \gamma _{10}(\phi ) \ne 0,\\ 0.1 &{} \text{ else }, \end{array}\right. } \\{} & {} \beta _m(\phi )= \displaystyle 4e^{\gamma _0(\phi )/1.8}, \quad \alpha _m(\phi ) = {\left\{ \begin{array}{ll} \displaystyle \gamma _{25}(\phi )\left( e^{\gamma _{25}(\phi )} - 1\right) ^{-1} &{} \text{ if } \; \displaystyle \gamma _{25}(\phi ) \ne 0,\\ 1 &{} \text{ else } , \end{array}\right. } \\{} & {} \beta _h(\phi ) = \left( \displaystyle e^{\gamma _{30}(\phi )} + 1\right) ^{-1}, \quad \alpha _h(\phi ) = 0.07e^{\gamma _0(\phi )/2}. \end{aligned}$$

We use the same rate functions as in Pods et al. (2013) (under the assumption that the temperature factor $C_T$ is equal to 1) which are themselves equivalent to the rate functions used in Lopreore et al. (2008). The notation is however modified to describe the resting potential $\phi _r$ as a parameter that could eventually vary and to give explicitly the value of the rate functions when the denominators are equal to zero. Computing the flux in (7) requires to solve the ODE system (11).