Reduced order models of myelinated axonal compartments

Abstract

The paper presents a hierarchical series of computational models for myelinated axonal compartments. Three classes of models are considered, either with distributed parameters (2.5D EQS–ElectroQuasi Static, 1D TL-Transmission Lines) or with lumped parameters (0D). They are systematically analyzed with both analytical and numerical approaches, the main goal being to identify the best procedure for order reduction of each case. An appropriate error estimator is proposed in order to assess the accuracy of the models. This is the foundation of a procedure able to find the simplest reduced model having an imposed precision. The most computationally efficient model from the three geometries proved to be the analytical 1D one, which is able to have accuracy less than 0.1%. By order reduction with vector fitting, a finite model is generated with a relative difference of 10^− 4 for order 5. The dynamical models thus extracted allow an efficient simulation of neurons and, consequently, of neuronal circuits. In such situations, the linear models of the myelinated compartments coupled with the dynamical, non-linear models of the Ranvier nodes, neuronal body (soma) and dendritic tree give global reduced models. In order to ease the simulation of large-scale neuronal systems, the sub-models at each level, including those of myelinated compartments should have the lowest possible order. The presented procedure is a first step in achieving simulations of neural systems with accuracy control.

Appendix: Derivation of the analytical (2.5D) model

The equation div(σ gradV ) = 0, satisfied by the potential V in each homogeneous subdomain, has the following form in cylindrical coordinates:

$$ \text{div}\left( \sigma\text{grad}V\right)=\frac{1}{r}\frac{\partial}{\partial r}\left( r\sigma\frac{\partial V}{\partial r}\right) + \frac{\partial}{\partial x}\left( \sigma\frac{\partial V}{\partial x}\right) = 0, $$

(26)

where the axial symmetry of the function (i.e. independence on the azimuthal angle) was taken into account. In each homogeneous subdomain the potential V is a harmonic function, a solution of the Laplace equation. According to the separation of variables method (Ioan 1988), in each homogeneous subdomain the potential V is assumed to have the form:

$$ V(x,r)=X(x)R(r). $$

(27)

Substituting Eq. (27) in Eq. (26) it follows that the partial differential equation can be decomposed into two ordinary linear differential equations satisfied by the two functions X and R:

$$ -\frac{1}{r}\frac{(\sigma rR^{\prime})^{\prime}}{\sigma R} = \frac{X^{\prime\prime}}{X} \Rightarrow \left\{\begin{array}{ll} \frac{X^{\prime\prime}}{X} = \lambda^{2} \\ -\frac{(\sigma rR^{\prime})'}{\sigma rR} = \lambda^{2}, \end{array}\right. $$

(28)

where λ is a positive real constant, called constant of separation.

The solution of the first equation is:

$$ X(x) = A \text{sh}(\lambda x) + B \text{ch}(\lambda x), $$

(29)

where A and B are integration constants.

In Eq. (28) σ is piecewise constant, (σ₁ for 0 < r < a and σ₂ for a < r < b). Therefore, on each homogeneous subdomain, the function R is the solution of the differential equation:

$$ r^{2}R^{\prime\prime} + rR^{\prime} + \lambda^{2}r^{2}R = 0, $$

which is a combination of zero order Bessel functions, with the following form for the general solution:

$$ R(r) = CJ_{0}(\lambda x) + DY_{0}(\lambda x). $$

The parameter λ and the integration constants are derived by imposing the boundary conditions:

• on terminal 2 (x = L, 0 < r < a):

  $$ \left. \frac{\text{d} V}{\text{d} n}\right|_{x=L} = 0 \Rightarrow X(x)=B^{\prime}\text{ch}(\lambda(L-x)), $$

  V (0,x) needs to have a finite value, so in the first subdomain (0 < r < a) D has to be zero and thus:

  $$ V(r,x) = CJ_{0}(\lambda x)\text{ch}(\lambda(L-x)). $$

  (30)

  After renaming the constants above, the potential has the general form:

  $$ V(r,x)=\left\{\begin{array}{ll} BJ_{0}(\lambda r)\text{ch}(\lambda(L-x)), & 0{<}r{<}a, \\ (CJ_{0}(\lambda r) + DY_{0}(\lambda r))\text{ch}(\lambda(L-x)), & a{<}r{<}b. \end{array}\right. $$

• on the interface r = a:

  $$ \begin{array}{@{}rcl@{}} &&V_{1}(a,x)=V_{2}(a,x), 0<x<L \\ &&\left. \frac{\partial\underline{V_{1}}}{\partial n}\right|_{r=a} = \beta\left. \frac{\partial \underline{V_{2}}}{\partial n}\right|_{r=a}, \beta=\frac{\sigma_{2}+j\omega\varepsilon_{2}}{\sigma_{1}+j\omega\varepsilon_{1}} = \frac{\sigma_{2}}{\sigma_{1}}, \text{for } \omega=0. \end{array} $$

  It follows that:

  $$ \left\{\begin{array}{ll} \frac{D}{C}=\left( \frac{1-\beta}{\beta}\right)/\left( \frac{Y_{1}(\lambda a)}{J_{1}(\lambda a)} - \frac{1}{\beta}\frac{Y_{0}(\lambda a)}{J_{0}(\lambda a)}\right); \\ \frac{B}{C}=1+\frac{D}{C}\frac{Y_{0}(\lambda a)}{J_{0}(\lambda a)} = 1+\frac{(1-\beta)J_{1}(\lambda a)Y_{0}(\lambda a)}{\beta J_{0}(\lambda a)Y_{1}(\lambda a) - J_{1}(\lambda a)Y_{0}(\lambda a)}. \end{array}\right. $$

• on the boundary r = b:

$$ V(b,x)=0, \forall x {\in} [0,L] \Rightarrow CJ_{0}(\lambda b) + DY_{0}(\lambda b) = 0. $$

This leads to the eigenvalues equation:

$$ \begin{array}{@{}rcl@{}} &&(1-\beta)Y_{0}(\lambda b)J_{0}(\lambda a)J_{1}(\lambda a) + \\ &&J_{0}(\lambda b)(\beta Y_{1}(\lambda a)J_{0}(\lambda a) - Y_{0}(\lambda a)J_{1}(\lambda a) ) = 0, \end{array} $$

which has an infinite number of solutions λ_k,k = 1, 2,..., $\infty$ and the general solution of the problem V (r,x) is obtained by superposition of all possible general forms:

$$ V(r,x){=}\left\{\begin{array}{ll} {\sum}_{k}{C_{k}\frac{B_{k}}{C_{k}}J_{0}(\lambda_{k}{r})\text{ch}(\lambda_k(L-x))}, 0<r<a \\ {\sum}_{k}{C_{k}(J_{0}(\lambda_{k}{r}) {+} \frac{D_{k}}{C_{k}}Y_{0}(\lambda_{k}{r}))\text{ch}(\lambda_k(L-x))}. \end{array}\right. $$

Equivalently:

$$ V(r,x)={\sum}_{k}{C_{k}R(\lambda_{kr})\text{ch}(\lambda_k(L-x))} , $$

where R are eigenfunctions given by:

$$ R(\lambda_{kr})=\left\{\begin{array}{ll} \frac{B_{k}}{C_{k}}J_{0}(\lambda_{k}{r}), & 0<r<a \\ J_{0}(\lambda_{k}{r}) + \frac{D_{k}}{C_{k}}Y_{0}(\lambda_{k}{r}). & a<r<b \end{array}\right. $$

The constant C_k is computed by imposing the Neumann boundary condition at x = 0:

$$ \begin{array}{@{}rcl@{}} \left. \frac{\partial V}{\partial x}\right|_{x=0} &=& f(r), \\ f(r) &=& -\frac{I_{1}}{\sigma\pi a^{2}}h(a-r)=\left\{\begin{array}{ll} -\frac{I_{1}}{\sigma\pi a^{2}}, & r \in (0,a) \\ 0, & r \in (a,b), \end{array}\right. \end{array} $$

(31)

where h is the Heaviside function (unit step). The function f(r) can be expanded into Fourier-Bessel series of eigenfunctions:

$$ f(r){=}{\sum}_{k=1}^{\infty}{F_{k}R(\lambda_{kr})}{=}{\sum}_{k=1}^{\infty}{C_{k}R(\lambda_{kr})\lambda_{k}\text{sh}(\lambda_{k} L)}, $$

(32)

where the Fourier coefficients F_k of this series result from the orthogonality property of the eigenfunctions:

$$ {<}R(\lambda_{j} r),R(\lambda_{k} r){>}{=}{{\int}_{0}^{b}}{r\sigma(r)R(\lambda_{j} r)R(\lambda_{k} r)}\text{d} r = 0, j \ne k $$

In EC regime and with j = k this relation becomes:

$$ \begin{array}{@{}rcl@{}} \|R_{k}\|^{2} &=& \sigma_{1} {{\int}_{0}^{a}}{r\frac{{B^{2}_{k}}}{{C^{2}_{k}}}J^{2}_{0}(\lambda_{k}{r})}\text{d} r + \\ &&\sigma_{2} {{\int}_{a}^{b}}{r(J_{0}(\lambda_{k}{r}) + \frac{D_{k}}{C_{k}}Y_{0}(\lambda_{k}{r}))^{2}}\text{d} r \Rightarrow \end{array} $$

$$ \begin{array}{@{}rcl@{}} \|R_{k}\|^{2} &=& \sigma_{1} \frac{{B^{2}_{k}}}{{C^{2}_{k}}} \frac{a^{2}}{2} (J^{2}_{0}(\lambda_{k}{a}) + J^{2}_{1}(\lambda_{k}{a})) \\ &&+ \sigma_{2} \frac{1}{2} (b^{2}(J^{2}_{0}(\lambda_{k}{b}) + J^{2}_{1}(\lambda_{k}{b})) \\ &&- a^{2}(J^{2}_{0}(\lambda_{k}{a}) + J^{2}_{1}(\lambda_{k}{a}))) \\ &&+ \sigma_{2} \frac{{D^{2}_{k}}}{{C^{2}_{k}}} \frac{1}{2} (b^{2}(Y^{2}_{0}(\lambda_{k}{b}) + Y^{2}_{1}(\lambda_{k}{b})) \\&&- a^{2}(Y^{2}_{0}(\lambda_{k}{a}) + Y^{2}_{1}(\lambda_{k}{a})) ) \\ &&+ \sigma_{2} \frac{D_{k}}{C_{k}} (b^{2}(J_{0}(\lambda_{k}{b})Y_{0}(\lambda_{k}{b}) \\ &&+ J_{1}(\lambda_{k}{b})Y_{0}(\lambda_{k}{b})) - a^{2}(J_{0}(\lambda_{k}{a})Y_{0}(\lambda_{k}{a}) \\ &&+ J_{1}(\lambda_{k}{a})Y_{0}(\lambda_{k}{a}) ) ). \end{array} $$

(33)

The scalar product between R(λ_kr) and f(r) leads to the expression of F_k:

$$ \begin{array}{@{}rcl@{}} &&<R(\lambda_{k} r),f(r)>=F_{k}\|R_{k}\|^{2} \Rightarrow \\ &&F_{k} {=} \frac{I_{1}}{\pi a^{2}\|R_{k}\|^{2}} {{\int}_{0}^{a}}rR(\lambda_{kr})\text{d} r {=} \frac{I_{1}}{\pi a\|R_{k}\|^{2}}\frac{B_{k}}{C_{k}}\frac{1}{\lambda_{k}}J_{1}(\lambda_{k}{a}).\\ \end{array} $$

(34)