Merge pull request #33 from halfflat/feature/cable-computation

Add analytic simple cable solver.

docs/passive_cable/cable.bib

+@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}
+}

docs/passive_cable/cable_computation.tex

+\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}

docs/passive_cable/makefile

+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

scripts/PassiveCable.jl

+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

scripts/README.md

-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.

scripts/profstats

-#!env python
+#!/usr/bin/env python2
 #coding: utf-8
 
 import json

scripts/tsplot

-#!env python
+#!/usr/bin/env python2
 #coding: utf-8
 
 import argparse