diff --git a/docs/passive_cable/cable.bib b/docs/passive_cable/cable.bib new file mode 100644 index 0000000000000000000000000000000000000000..1f47e75b738ede826221dbc4669e419c21bfaa30 --- /dev/null +++ b/docs/passive_cable/cable.bib @@ -0,0 +1,33 @@ +@article{bhalla1992, + title = {Rallpacks: a set of benchmarks for neuronal simulators}, + author = {Upinder S. Bhalla and David H. Bilitch and James M. Bower}, + journal = {Trends in Neurosciences}, + volume = {15}, + number = {11}, + pages = {453--458}, + year = {1992}, + month = {11}, + doi = {10.1016/0166-2236(92)90009-W} +} + +@article{churchill1937, + title = {The inversion of the Laplace transformation by a direct expansion in series and its application to boundary-value problems}, + author = {R. V. Churchill}, + journal = {Mathematische Zeitschrift}, + volume = {42}, + number = {1}, + pages = {567--579}, + year = {1937}, + doi = {10.1007/BF01160095} +} + +@article{borwein2009, + title = {Uniform bounds for the complementary incomplete gamma function}, + author = {Jonathan M Borwein and O-Yeat Chan}, + journal = {Mathematical Inequalities and Applications}, + volume = {12}, + number = {1}, + pages = {115--121}, + year = {2009}, + doi = {10.7153/mia-12-10} +} diff --git a/docs/passive_cable/cable_computation.tex b/docs/passive_cable/cable_computation.tex new file mode 100644 index 0000000000000000000000000000000000000000..39a372196c006c7c67837f039df88c903cb439b4 --- /dev/null +++ b/docs/passive_cable/cable_computation.tex @@ -0,0 +1,324 @@ +\documentclass[parskip=half]{scrartcl} + +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage[citestyle=authoryear-icomp,bibstyle=authoryear]{biblatex} +\usepackage{booktabs} +\usepackage[style=british]{csquotes} +\usepackage{fontspec} +\usepackage{siunitx} + +\setmainfont[Numbers=Lowercase]{TeX Gyre Pagella} +\newfontfamily\dispfamily{TeX Gyre Pagella} +\addtokomafont{disposition}{\dispfamily} + +\MakeAutoQuote{«}{»} +\DeclareMathOperator{\Res}{Res} + +\addbibresource{cable.bib} + +\title{Analytic solution to passive cable model} +\author{Sam Yates} +\date{October 20, 2016} + +\begin{document} + +\maketitle + +\section{Background} + +The Rallpack suite \parencite{bhalla1992} is a set of three tests for cable-based neuronal +simulators, for validation and benchmarking. Validation data is provided with the +distribution of the suite, which at the time of writing can now be found within the +\textsc{Genesis} simulator simulator distribution.\footnote{see \url{http://genesis-sim.org}.} + +The first two of the tests are based on passive cable models which admit analytic +solutions, and the Rallpack suite includes code that can generate the reference curves +for these models. The reference code, however, requires hand-tuning of the iterations; +the authors state in the internal documentation that +«the use of 10 terms gives at least six figure accuracy», but this is for the +fixed time increments used in their dataset and is not validated within the code +or by the supplied material. + +For flexibility of testing and cross-checking of the Rallpack reference data, +it would be useful to have at hand a robust calculator of the passive cable +potential solution. + +\section{The Rallpack 1 model} + +The model comprises a cylindrical cable segment with the physical and electrical +properties as described in Table~\ref{tbl:rallpack1}. A current of \SI{0.1}{\nA} +is injected at the left end of the cable from $t=0$, and the initial potential +is the reversal potential. + +\begin{table}[ht] + \centering + \begin{tabular}{lSl} + \toprule + Term & {Value} & Property\\ + \midrule + $d$ & \SI{1.0}{\um} & cable diameter \\ + $L$ & \SI{1.0}{\mm} & cable length \\ + $R_A$ & \SI{1.0}{\ohm\m} & bulk axial resistivity \\ + $R_M$ & \SI{4.0}{\ohm\m\squared} & areal membrane resistivity \\ + $C_M$ & \SI{0.01}{\F\per\m\squared} & areal membrane capacitance \\ + $E_M$ & \SI{-65.0}{\mV} & membrane reversal potential \\ + \bottomrule + \end{tabular} + \caption{Cable properties for the Rallpack 1 model.} + \label{tbl:rallpack1} +\end{table} + +The potential on a constant-radius passive cable is governed by the PDE +\begin{equation} + \lambda^2 \frac{\partial^2 v}{\partial x^2} = + \tau\frac{\partial v}{\partial t} + v - E, +\end{equation} +where $E$ is the membrane reversal potential, and $\lambda$ and $\tau$ are the +electrotonic length and time constants for the cable. It is convenient +to consider the electrical properties of the cable per unit-length, in terms +of the linear axial resistivity $r$, the linear membrane capacitance $r$, +and the linear membrane capacitance $c$. These determine $\lambda$ and $\tau$ by +\begin{gather*} + \lambda = \frac{1}{r\cdot g},\\ + \tau = r\cdot c. +\end{gather*} + +With the model boundary conditions, +\begin{subequations} + \begin{align} + v(x, 0) &= E, \\ + \left.\frac{\partial v}{\partial x}\right\vert_{x=0} & = -Ir, \\ + \left.\frac{\partial v}{\partial x}\right\vert_{x=L} & = 0, + \end{align} +\end{subequations} +where $I$ is the injected current and $L$ is the cable length. + +The solution $v(x, t)$ can be expressed in terms of the solution $g(x, t; L)$ +to a normalized version of the cable equation, +\begin{subequations} + \begin{align} + \label{eq:normcable} + \frac{\partial^2 g}{\partial x^2} & = + \frac{\partial g}{\partial t} + g, + \\ + \label{eq:normcableinitial} + g(x, 0) &= 0, + \\ + \label{eq:normcableleft} + \left.\frac{\partial g}{\partial x}\right\vert_{x=0} & = 1, + \\ + \label{eq:normcableright} + \left.\frac{\partial g}{\partial x}\right\vert_{x=L} & = 0 + \end{align} +\end{subequations} +by +\begin{equation} + v(x, t)= E - Ir\lambda \cdot g(\frac{x}{\lambda}, \frac{t}{\tau}; \frac{L}{\lambda}). +\end{equation} + +\section{Solution to the normalized cable equation} + +Let $G(x, s)$ be the Laplace transform of $g(x, t)$ with respect to $t$. +If the inverse transform of $G$ exists, it must agree with $g$ almost everywhere +for $t>0$, or for all $t>0$ given the continuity of $g$. + +From +\eqref{eq:normcable} and \eqref{eq:normcableinitial}, +\begin{equation} + \label{eq:lap} + \frac{\partial^2 G}{\partial x^2} = sG + G. +\end{equation} +The boundary conditions \eqref{eq:normcableleft} and \eqref{eq:normcableright} give +\begin{align} + \label{eq:lapleft} + \frac{\partial G}{\partial x}(0, s) & = \frac{1}{s} + \\ + \label{eq:lapright} + \frac{\partial G}{\partial x}(L, s) & = 0. +\end{align} + +Solutions to \eqref{eq:lap} must be of the form +\begin{equation} + G(x, s) = A(s) e^{mx} +B(s) e^{-mx} +\end{equation} +where $m=\sqrt{1+s}$. From \eqref{eq:lapleft} and \eqref{eq:lapright}, +\begin{align*} + mA(s) - mB(s) & = \frac{1}{s} + \\ + mA(s)e^{mL} -mB(s)e^{-mL} & = 0 +\end{align*} +and thus +\begin{align*} + A(s) & = \frac{1}{sm (1-e^{2mL})} + \\ + B(s) & = \frac{e^{2mL}}{sm (1-e^{2mL})}. +\end{align*} + +Consequently, +\begin{equation} + \begin{aligned} + G(x, s) &= \frac{1}{ms}\cdot\frac{e^{mx}+e^{2mL-mx}}{1-e^{2mL}}\\ + &= - \frac{1}{ms}\cdot\frac{\cosh m(L-x)}{\sinh mL}. + \end{aligned} +\end{equation} + +\subsection*{Inversion} + +Sufficient conditions for the inverse transform of $G(x,s)$ to exist +and be representable in series form are as follows +\parencite[][Theorem 4]{churchill1937}: +\begin{enumerate} +\item + $G(x, s)$ is analytic in some right half-plane $H$, +\item + the singularities of $G(x, s)$ are all poles, and +\item for some $k>1$, $|s^k G(x, s)|$ is bounded in $H$ and + on the circles $|s|=\rho_i$, where $\rho_i$ is an unbounded monotonically + increasing sequence. +\end{enumerate} +If in addition all the poles $\sigma_j$ of $G$ are simple, then +the series expansion of the inverse transform is given by +\begin{equation} + g(x,t) = \sum_j e^{\sigma_j t}\, \Res(G;\sigma_j),\quad t>0. +\end{equation} + +$G(x,s)$ has poles when $s=0$, $m=0$, or $\sinh mL=0$. +In the following, let +\begin{gather*} + a_k = k\pi /L\\ + m_k = a_k i\\ + s_k = m_k^2 -1 = -k^2\pi^2/L^2-1. +\end{gather*} +for $k\in\mathbb{Z}$. +The poles are then $0$, $-1$, and $s_k$ for $k\geq 1$. + +There are no branch points +arising from $m=\sqrt{1+s}$, as letting $m=-\sqrt{1+s}$ leaves $G$ unchanged. + +For $|s+1|>\epsilon$, +\begin{equation} + \begin{aligned} + |s^{3/2}G(x,s)|^2 + & \leq (1+\epsilon)^{-1} \left| \frac{\cosh m(L-x)}{\sinh mL} \right|^2 + \\ + & \leq (1+\epsilon)^{-1} (1+|\coth mL|)^2 + \label{eq:gbounds} + \end{aligned} +\end{equation} +which is bounded in the half-plane $\Re(s) \geq -1+\epsilon$. +As $\coth u+iv$ is periodic in $v$, \eqref{eq:gbounds} is also bounded outside +the regions $\{s: m\in D_k\}$, where $D_k$ are non-overlapping $\delta$-neighborhoods of $m_k$ +for some fixed small $\delta>0$. These are contained within disks of radius +$2k\pi\delta/L+\delta^2$ about $s_k$ for $k\geq 1$, which then are separated by +circles $|s|=\rho_k$ about the origin for some $\rho_k$ in $(|s_k|, |s_k+1|)$. + +\textbf{Pole at $s=0$}.\\ +Recalling $m=\sqrt{1+s}$, $m$ and $\sinh mL$ are non-zero in +a neighbourhood of $s=0$, and so the pole is simple and +\begin{equation} + \begin{aligned} + \Res(G; 0) & = - \frac{1}{m}\cdot\left.\frac{\cosh m(L-x)}{\sinh mL}\right|_{s=0}\\ + & = - \frac{\cosh (L-x)}{\sinh L}. + \end{aligned} +\end{equation} + +\textbf{Pole at $s=-1$}.\\ +Let $G(x,s)=f(x,s)/h(s)$, where +\begin{equation} + \begin{aligned} + f(x, s) &= -\frac{1}{s}\cosh m(L-x)\\ + h(s) &= m\sinh mL. + \end{aligned} + \label{eq:hf} +\end{equation} +Noting that $dm/ds = \frac{1}{2}m^{-1}$, +\begin{equation} + \begin{aligned} + h'(s) &= \frac{1}{2}m^{-1}\sinh mL + \frac{1}{2}L\cosh mL \\ + &= \frac{1}{2}L + \frac{1}{2}L + O(m^2) \quad(m\to 0) \\ + &= L + O(s+1) \quad(s\to -1). + \label{eq:hprime} + \end{aligned} +\end{equation} +The pole is therefore simple, and +\begin{equation} + \Res(G; -1) = f(x, -1)/h'(-1) = -\frac{1}{L}. +\end{equation} + +\textbf{Pole at $s=s_k$, $k \geq 1$}.\\ +Taking $h$ and $f$ as above \eqref{eq:hf}, +$m_k$ is non-zero for $k\geq 1$ and +\[ + h'(s_k) = \frac{1}{2}L\cosh m_kL. +\] +Consequently the pole is simple and +\begin{equation} + \begin{aligned} + \Res(G; s_k) + & = f(x, s_k)/h'(s_k)\\ + & = -\frac{2}{s_k L}\frac{\cosh m_k(L-x)}{\cosh m_kL} \\ + & = -\frac{2}{s_k L}\frac{\cosh m_kL\cosh m_kx-\sinh m_kL\sinh m_kL}{\cosh m_kL} \\ + & = -\frac{2}{s_k L}\cosh m_k x, + \end{aligned} +\end{equation} +as $\sinh m_k=0$. + +In terms of $a_k$, +\begin{equation} + \Res(G; s_k) = \frac{2}{L}\cdot\frac{1}{1+a_k^2}\cdot\cos a_kx. +\end{equation} + +The series exapnsion for $g(x, t)$ therefore is +\begin{equation} + g(x, t) = -\frac{\cosh(L-x)}{\sinh L} + \frac{1}{L}e^{-t}\left\{ + 1+2\sum_{k=1}^\infty \frac{e^{-ta_k^2}}{1+a_k^2}\cos a_k x\right\}. + \label{eq:theg} +\end{equation} + +\section{Computation of series} + +As $t$ approaches zero, the series \eqref{eq:theg} converges increasingly slowly. +For computational purposes, an estimate on the series residual allows the determination +of stopping criteria for a given tolerance. + +Let $g_n$ be the partial sum +\begin{equation} + g_n(x, t) = -\frac{\cosh(L-x)}{\sinh L} + \frac{1}{L}e^{-t}\left\{ + 1+2\sum_{k=1}^n \frac{e^{-ta_k^2}}{1+a_k^2}\cos a_k x\right\}. +\end{equation} +so that $g(x, t) =\lim_{n\to\infty} g_n(x,t)$. Let $\bar{g}_n = |g-g_n|$ be the +residual. The $a_k$ form an increasing sequence, so +\begin{equation} + \begin{aligned} + \bar{g}_n(x,t) + & \leq \frac{2}{L}e^{-t}\sum_{n+1}^\infty\frac{e^{-ta_k^2}}{1+a_k^2}\\ + & \leq \frac{2}{L}e^{-t}\int_{a_n}^\infty \frac{e^{-tu^2}}{1+u^2}\,du\\ + & < \frac{2}{L}e^{-t}\int_{a_n}^\infty \frac{e^{-tu^2}}{u^2}\,du. + \end{aligned} + \label{eq:gbar} +\end{equation} +Substituting $u=\sqrt{v/t}$ gives the identity +\begin{equation} + \int_{a}^\infty \frac{e^{-tu^2}}{u^2}\,du + = \frac{1}{2}\sqrt{t}\int_{a^2t}^\infty e^{-tv} v^{\frac{-3}{2}}\,dv + = \frac{1}{2}\sqrt{t}\,\Gamma(-\frac{1}{2},a^2 t), +\end{equation} +where $\Gamma(\alpha, z)$ is the upper incomplete gamma function. + +For real $\alpha<1$ and $z>0$, \textcite[][Theorem 2.3]{borwein2009} give the upper bound +\[ + \Gamma(\alpha, z) \leq z^{\alpha-1} e^{-z}. +\] +Substituting into \eqref{eq:gbar} gives +\begin{equation} + \begin{aligned} + \bar{g}_n(x,t) + & < \frac{1}{L}e^{-t}\sqrt{t}\,\Gamma(-\frac{1}{2},a_n^2 t) \\ + & \leq \frac{1}{L}e^{-t}\sqrt{t}\,(a_n^2 t)^{-\frac{3}{2}}e^{-a_n^2t} \\ + & = \frac{e^{-t(1+a_n^2)}}{L t a_n^3}. + \end{aligned} +\end{equation} + +\printbibliography +\end{document} diff --git a/docs/passive_cable/makefile b/docs/passive_cable/makefile new file mode 100644 index 0000000000000000000000000000000000000000..8756af0b9b8f33a1a145a07259d07a335a86de86 --- /dev/null +++ b/docs/passive_cable/makefile @@ -0,0 +1,14 @@ +SOURCES := cable_computation.tex + +LATEXMK := latexmk -e '$$clean_ext=q/bbl run.xml/' -xelatex -use-make -halt-on-error + +all: cable_computation.pdf + +cable_computation.pdf: cable_computation.tex cable.bib + $(LATEXMK) $< + +clean: + for s in $(SOURCES); do $(LATEXMK) -c "$$s"; done + +realclean: + for s in $(SOURCES); do $(LATEXMK) -C "$$s"; done diff --git a/scripts/PassiveCable.jl b/scripts/PassiveCable.jl new file mode 100644 index 0000000000000000000000000000000000000000..76450490ce5c4f265403edd35f2bedad4ad7c3bb --- /dev/null +++ b/scripts/PassiveCable.jl @@ -0,0 +1,141 @@ +module PassiveCable + +export cable_normalize, cable, rallpack1 + +# Compute solution g(x, t) to +# +# ∂²g/∂x² - g - ∂g/∂t = 0 +# +# on [0, L] × [0,∞), subject to: +# +# g(x, 0) = 0 +# ∂g/∂x (0, t) = 1 +# ∂g/∂x (L, t) = 0 +# +# Parameters: +# x, t, L: as described above +# tol: absolute error tolerance in result +# +# Return: +# g(x, t) + +function cable_normalized(x, t, L; tol=1e-8) + if t<=0 + return 0.0 + else + ginf = -cosh(L-x)/sinh(L) + sum = exp(-t/L) + + for k = countfrom(1) + a = k*pi/L + e = exp(-t*(1+a^2)) + + sum += 2/L*e*cos(a*x)/(1+a^2) + resid_ub = e/(L*a^3*t) + + if resid_ub<tol + break + end + end + return ginf+sum + end +end + + +# Compute solution f(x, t) to +# +# λ²∂²f/∂x² - f - τ∂f/∂t = 0 +# +# on [0, L] x [0, ∞), subject to: +# +# f(x, 0) = V +# ∂f/∂x (0, t) = I·r +# ∂f/∂x (L, t) = 0 +# +# where: +# +# λ² = 1/(r·g) length constant +# τ = r·c time constant +# +# In the physical model, the parameters correspond to the following: +# +# L: length of cable +# r: linear axial resistivity +# g: linear membrane conductivity +# c: linear membrane capacitance +# V: membrane reversal potential +# I: injected axial current on the left end (x = 0) of the cable. +# +# Note that r, g and c are specific 1-d quantities that differ from +# the cable resistivity r_L, specific membrane conductivity ḡ and +# specific membrane capacitance c_m as used elsewhere. If the +# cross-sectional area is A and cable circumference is f, then +# these quantities are related by: +# +# r = r_L/A +# g = ḡ·f +# c = c_m·f +# +# Parameters: +# x: displacement along cable +# t: time +# L, lambda, tau, r, V, I: as described above +# tol: absolute error tolerance in result +# +# Return: +# computed potential at (x,t) on cable. + +function cable(x, t, L, lambda, tau, r, V, I; tol=1e-8) + scale = I*r*lambda; + if scale == 0 + return V + else + tol_n = abs(tol/scale) + return scale*cable_normalized(x/lambda, t/tau, L/lambda, tol=tol_n) + V + end +end + + +# Rallpack 1 test +# +# One sided cable model with the following parameters: +# +# RA = 1 Ω·m bulk axial resistivity +# RM = 4 Ω·m² areal membrane resistivity +# CM = 0.01 F/m² areal membrane capacitance +# d = 1 µm cable diameter +# EM = -65 mV reversal potential +# I = 0.1 nA injected current +# L = 1 mm cable length. +# +# (This notation aligns with that used in the Rallpacks paper.) +# +# Note that the injected current as described in the Rallpack model +# is trans-membrane, not axial. Consequently we need to swap the +# sign on I when passing to the cable function. +# +# Parameters: +# x: displacement along cable [m] +# t: time [s] +# tol: absolute tolerance for reported potential. +# +# Return: +# computed potential at (x,t) on cable. + +function rallpack1(x, t; tol=1e-8) + RA = 1 + RM = 4 + CM = 1e-2 + d = 1e-6 + EM = -65e-3 + I = 0.1e-9 + L = 1e-3 + + r = 4*RA/(pi*d*d) + lambda = sqrt(d/4 * RM/RA) + tau = CM*RM + + return cable(x, t, L, lambda, tau, r, EM, -I, tol=tol) +end + +end # module PassiveCable diff --git a/scripts/README.md b/scripts/README.md index e645a762eaf88c5f7196a995502589da44c2f89c..373f9f1330840c7cddf9fcf264d92f943c55c33b 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -1,5 +1,4 @@ -tsplot ------- +#tsplot The `tsplot` script is a wrapper around matplotlib for displaying a collection of time series plots. @@ -78,4 +77,58 @@ of the timeseries data. Use the `-o` or `--output` option to save the plot as an image, instead of displaying it interactively. +#profstats +`profstats` collects the profiling data output from multiple MPI ranks and performs +a simple statistical summary. + +Input files are in the JSON format emitted by the profiling code. + +By default, `profstats` reports the quartiles of the times reported for each +profiling region and subregion. With the `-r` option, the collated raw times +are reported instead. + +Output is in CSV format. + +#PassiveCable.jl + +Compute analytic solutions to the simple passive cylindrical dendrite cable +model with step current injection at one end from t=0. + +This is used to generate validation data for the first Rallpack test. + +Module exports the following functions: + + * `cable_normalized(x, t, L; tol)` + + Compute potential _V_ at position _x_ in [0, _L_] at time _t_ ≥ 0 according + to the normalized cable equation with unit length constant and time constant. + + Neumann boundary conditions: _V'_(0) = 1; _V'_(L) = 0. + Initial conditions: _V_( _x_, 0) = 0. + + Absolute tolerance `tol` defaults to 1e-8. + + * `cable(x, t, L, lambda, tau, r, V, I; tol)` + + Compute the potential given: + * length constant `lambda` + * time constant `tau`, + * axial linear resistivity `r` + * injected current of `I` at the origin + * reversal potential `V` + + Implied units must be compatible, e.g. SI units. + + Absolute tolerance `tol` defaults to 1e-8. + + * `rallpack1(x, t; tol)` + + Compute the value of the potential in the Rallpack 1 test model at position + `x` and time `t`. + + Parameters for the underlying cable equation calculation are taken from the + Rallpack model description in SI units; as the cable length is 1 mm in this + model, `x` can take values in [0, 0.001]. + + Absolute tolerance `tol` defaults to 1e-8. diff --git a/scripts/profstats b/scripts/profstats index 33f37356836104e57758f0407e1da5fd0625a111..88f68c72e6253fa4c0d240d051d7ce1be1960604 100755 --- a/scripts/profstats +++ b/scripts/profstats @@ -1,4 +1,4 @@ -#!env python +#!/usr/bin/env python2 #coding: utf-8 import json diff --git a/scripts/tsplot b/scripts/tsplot index 9ad6e2a5c714a5f634e37c31c88411f63923f130..8e89cbda96451680bfc9de8087ba558c0282defc 100755 --- a/scripts/tsplot +++ b/scripts/tsplot @@ -1,4 +1,4 @@ -#!env python +#!/usr/bin/env python2 #coding: utf-8 import argparse