From c2a84a005fd01ad7bb1106e2ef6c9d8bcf3b7ccd Mon Sep 17 00:00:00 2001
From: bcumming <bcumming@cscs.ch>
Date: Tue, 19 Jan 2016 17:17:07 +0100
Subject: [PATCH] refine the FVM formulation
---
docs/formulation.tex | 53 +++++++++++++++++++++++--------------------
docs/images/cable.tex | 4 ++--
docs/report.tex | 3 ++-
docs/symbols.tex | 24 ++++++++++++++++++++
4 files changed, 57 insertions(+), 27 deletions(-)
diff --git a/docs/formulation.tex b/docs/formulation.tex
index 755817f7..bf8f2fe4 100644
--- a/docs/formulation.tex
+++ b/docs/formulation.tex
@@ -48,7 +48,7 @@ where
\item $c_m$ is the {specific membrane capacitance}, approximately the same for all neurons $\approx 10~nF/mm^2$. Related to \emph{membrane capacitance} $C_m$ by the relationship $C_m=c_{m}A$, where $A$ is the surface area of the cell.
\item $i_m$ is the membrane current $[A]$. The total contribution from ion and synaptic channels is expressed as a the product of current per unit area $i_m$ and the surface area.
\item $i_e$ is the electrode current flowing into the cell, divided by surface area, i.e. $i_e=I_e/A$.
- \item $r_L$ is intracellular resistivity, typical value $1~k\Omega$
+ \item $r_L$ is intracellular resistivity, typical value $1~k\Omega \text{cm}$
\end{itemize}
Note that the standard convention is followed, whereby membrane and synapse currents ($i_m$) are positive when outward, and electrod currents ($i_e$) are positive inward.
@@ -62,8 +62,8 @@ The PDE in (\ref{eq:cable}) is derived from the following mass balance expressio
%\end{align}
\begin{equation}
\int_{\Omega_i}{c_m \pder{V}{t} } \deriv{v} =
- + \sum_{j\in\mathcal{N}_i} {\int_{\Gamma_{i,j}} q_{i,j} \deriv{s} }
- + \int_{\Gamma_{ext}} {q_i} \deriv{s}
+ - \sum_{j\in\mathcal{N}_i} {\int_{\Gamma_{i,j}} q_{i,j} \deriv{s} }
+ - \int_{\Gamma_{i}} {q_i} \deriv{s}
\label{eq:cable_balance}
\end{equation}
where
@@ -71,12 +71,12 @@ where
\item $\int_\Omega \cdot \deriv{v}$ is shorthand for the volume integral over the segment $\Omega_i$
\item $\int_\Gamma \cdot \deriv{s}$ is shorthand for the surface integral over the surface $\Gamma$
\item $q_{i,j}=-\frac{1}{r_L}\pder{V}{x} n_{i,j}$ is the flux per unit area of current \emph{from segment $i$ to segment $j$} over the interface $\Gamma_{i,j}$ between the two segments.
- \item $q_i=i_m - i_e$ is the flux per unit area over the cell membrane $\Gamma_i$ due to ion channels, synapses and electrodes.
+ \item $q_i=i_m - i_e$ is the flux per unit area over the cell membrane $\Gamma_i$ due to ion channels, synapses and electrodes (where $q_i>0$ implies flux out of the cell).
\item the set $\mathcal{N}_i$ is the set of segments that are neighbours of $\Omega_i$
\end{itemize}
The surface of the cable segment is sub-divided into the internal and external surfaces.
-The external surface $\Gamma_{ext}$ is the cell membrane at the interface between the extra-cellular and intra-cellular regions.
+The external surface $\Gamma_{i}$ is the cell membrane at the interface between the extra-cellular and intra-cellular regions.
The current, which is the conserved quantity in our conservation law, over the surface is composed of the synapse and ion channel contributions.
This is derived from a thin film approximation to the cell membrane, whereby the membrane is treated as an infinitesimally thin interface between the intra and extra cellular regions.
@@ -94,12 +94,12 @@ If the cell is represented by cylinders or frustrums\footnote{a frustrum is a tr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Finite volume discretization}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-The finite volume method is a natural choice for the solution of the conservation law in~\eq{eq:cable_balance} and~\eq{eq:cable_balance_1D}.
+The finite volume method is a natural choice for the solution of the conservation law in~\eq{eq:cable_balance}.
\begin{itemize}
\item the $x_i$ are spaced uniformly with distance $x_{i+1}-x_{i} = \Delta x$
\item control volumes are formed by locating the boundaries between adjacent points at $(x_{i+1}+x_{i})/2$
- \item this discretization differs from the finite difference method used in Neuron because the equation is explicitly solved for at the end of cable segments, and because the finite volume discretization is applied to all points. Neuron uses special algebraic \emph{zero area} formulation for nodes at branch points.
+ \item this discretization differs from the finite differences used in Neuron because the equation is explicitly solved for at the end of cable segments, and because the finite volume discretization is applied to all points. Neuron uses special algebraic \emph{zero area} formulation for nodes at branch points.
\end{itemize}
%-------------------------------------------------------------------------------
@@ -109,10 +109,10 @@ We proceed by defining the \emph{volume average} of a quantity $\varphi$ as foll
\begin{equation}
\bar{\varphi}_i = \frac{1}{\Delta_i} \int_{\Omega_i}{\varphi}\deriv{v},
\end{equation}
-where $\Omega_i$ is the control volume with the point $x_i$ at its centroid illustrated in \fig{fig:segment}, and $v_i$ is the volume of the control volume $\Omega_i$.
+where $\Omega_i$ is the control volume with the point $x_i$ at its centroid illustrated in \fig{fig:segment}, and $\Delta_i$ is the volume of the control volume $\Omega_i$.
The integral one the left hand side of~\eq{eq:cable_balance} can be expressed in terms of the volume average of $V$
\begin{equation}
- \int_{\Omega}{c_m \pder{V}{t} } \deriv{v} = \frac{c_m}{\Delta_i} \pder{\bar{V}_i}{t}.
+ \int_{\Omega}{c_m \pder{V}{t} } \deriv{v} = \Delta_i c_m \pder{\bar{V}_i}{t}.
\label{eq:dvdt}
\end{equation}
@@ -122,22 +122,25 @@ This is effectively treats voltage as a piecewise continuous funtion, with disco
%-------------------------------------------------------------------------------
\subsubsection{Intra-cellular flux}
%-------------------------------------------------------------------------------
-The intracellular flux terms in~\eq{eq:cable_balance} are sum of the flux over interface between adjacent cable segments
+The intracellular flux terms in~\eq{eq:cable_balance} are sum of the flux over the interfaces between compartment $i$ and its set of neighbouring compartments $\mathcal{N}_i$ is
\begin{equation}
- \sum_{j\in\mathcal{N}_i} {\int_{\Gamma_{i,j}} \left( \frac{1}{r_L}\pder{V}{x} n_{i,j} \right) \deriv{s} }
+ \sum_{j\in\mathcal{N}_i} { \int_{\Gamma_{i,j}} { q_{i,j} \deriv{s} } }.
\end{equation}
-
-For the segment centred at $x_i$ with neighbour $x_j$
+where the flux per unit area from compartment $i$ to compartment $j$ is
\begin{align}
- q_{i,j} &= \int_{\Gamma_{i,j}} \left( \frac{1}{r_L}\pder{V}{x} n_{i,j} \right) \deriv{s} \nonumber \\
- &= \int_{\Gamma_{i,j}} \frac{1}{r_L}\frac{V_i-V_j}{\Delta x_{i,j}} \deriv{s}
- \label{eq:q_ij_intermediate}
+ q_{i,j} = - \frac{1}{r_L}\pder{V}{x} n_{i,j}.
+ \label{eq:q_ij}
+\end{align}
+The derivative with respect to the outward-facing normal can be approximated as follows
+\begin{equation*}
+ \pder{V}{x} n_{i,j} \approx \frac{V_j - V_i}{\Delta x_{i,j}}
+\end{equation*}
+where $\Delta x_{i,j}$ is the distance between $x_i$ and $x_j$, i.e. $\Delta x_{i,j}=|x_i-x_j|$.
+Using this approximation for the derivative, the flux over the surface in~\eq{eq:q_ij} is approximated as
+\begin{align}
+ q_{i,j} \approx \frac{1}{r_L}\frac{V_i - V_j}{\Delta x_{i,j}}.
+ \label{eq:q_ij_intermediate}
\end{align}
-where:
-\begin{itemize}
- \item $q_{i,j}$ is the \emph{flux} of charge over the interface $\Gamma_{i,j}$
- \item $\Delta x_{i,j}$ is the distance between $x_i$ and $x_j$, i.e. $\Delta x_{i,j}=|x_i-x_j|$.
-\end{itemize}
The terms inside the integral in equation~\eq{eq:q_ij_intermediate} are constant everywhere on the surface $\Gamma_{i,j}$, so the integral becomes
\begin{align}
@@ -147,7 +150,7 @@ The terms inside the integral in equation~\eq{eq:q_ij_intermediate} are constant
&= \frac{\pi a_{i,j}^2}{r_L \Delta x_{i,j}} (V_i-V_j)
\label{eq:q_ij}
\end{align}
-where $\sigma_{i,j}=\pi a_{i,j}^2$ is the area of the surface $\Gamma_{i,j}$, which is a circle of radius $a_{r,j}$.
+where $\sigma_{i,j}=\pi a_{i,j}^2$ is the area of the surface $\Gamma_{i,j}$, which is a circle of radius $a_{i,j}$.
Some symmetries
\begin{itemize}
@@ -187,9 +190,10 @@ where $a_{i,\ell}$ and $a_{i,r}$ are the radii of at the left and right end of t
\subsubsection{Putting it all together}
%-------------------------------------------------------------------------------
By substituting the volume averaging of the temporal derivative in~\eq{eq:dvdt} approximations for the flux over the surfaces in~\eq{eq:q_ij} and~\eq{eq:cv_volume} respectively into the conservation equation~\eq{eq:cable_balance} we get the following ODE defined for each node in the cell
+\note{there is an error somewhere, because the units on the lhs and the first summation term on the rhs do not match: $[lhs]=F\cdot V\cdot s^{-1}\cdot cm$ and $[rhs]=F\cdot V\cdot s^{-1}$ }
\begin{equation}
- \frac{c_m}{\Delta_i} \dder{\bar{V}_i}{t}
- = -\sum_{j\in\mathcal{N}_i} {\frac{\sigma_{i,j}}{r_L \Delta x_{i,j}} (V_i-V_j)} + \sigma_i\cdot(i_m(V_i) - i_e(x_i)),
+ \Delta_i c_m \dder{V_i}{t}
+ = -\sum_{j\in\mathcal{N}_i} {\frac{\sigma_{i,j}}{r_L \Delta x_{i,j}} (V_i-V_j)} - \sigma_i\cdot(i_m(V_i) - i_e(x_i)),
\label{eq:ode}
\end{equation}
where
@@ -226,3 +230,4 @@ The value of the lateral area and volume, $\sigma_i$ and $\Delta_i$ in~\eq{eq:si
\begin{equation}
\sigma_i = \sum_{j\in\mathcal{N}_i} {}
\end{equation}
+
diff --git a/docs/images/cable.tex b/docs/images/cable.tex
index f9008d0f..c11e2b74 100644
--- a/docs/images/cable.tex
+++ b/docs/images/cable.tex
@@ -70,8 +70,8 @@
\draw [] (Tface) -- ( -1, 1);
\draw [] (Bface) -- ( -1, 1);
\draw [] (1, -0.5) -- (1.5, -1);
-\node [] at (-2.5, 1.15) {$\Gamma_{{left}}$};
-\node [] at ( 2.5, 1.15) {$\Gamma_{{right}}$};
+\node [] at (-2.5, 1.15) {$\Gamma_{{i,i-1}}$};
+\node [] at ( 2.5, 1.15) {$\Gamma_{{i,i+1}}$};
\node [] at ( -1, 1.15) {$\Gamma_{{e}}$};
\node [] at ( 1.7, -1.2) {$\Omega$};
diff --git a/docs/report.tex b/docs/report.tex
index 24d7f765..54d93c8d 100644
--- a/docs/report.tex
+++ b/docs/report.tex
@@ -60,6 +60,7 @@
%----------------------------------------------------------------------------------------
\newcommand{\todo}[1]{\textbf{\textcolor{blue}{TODO: #1}}} % add a comment to the article
+\newcommand{\note}[1]{\textbf{\textcolor{red}{NOTE: #1}}} % add a note to the article
\newcommand{\tbl}[1]{\textbf{Table \ref{#1}}\xspace}
\newcommand{\fig}[1]{\textbf{Figure \ref{#1}}\xspace}
@@ -94,7 +95,7 @@
\section{Formulation}
\input{formulation.tex}
-\section{Apendix}
+\section{Appendix}
\input{appendix.tex}
\section{Symbols and Units}
diff --git a/docs/symbols.tex b/docs/symbols.tex
index 814b7d4a..b97daf69 100644
--- a/docs/symbols.tex
+++ b/docs/symbols.tex
@@ -41,3 +41,27 @@
\end{center}
\caption{Symbols and quantities.}
\end{table*}
+%-------------------------------------------------------------------------------
+\subsubsection{Units}
+%-------------------------------------------------------------------------------
+Reffering to the cable equation first defined in~\eq{eq:cable}
+\begin{equation*}
+ c_m \pder{V}{t} = \frac{1}{2\pi a r_{L}} \pder{}{x} \left( a^2 \pder{V}{x} \right) - i_m + i_e,
+\end{equation*}
+
+If the units are taken to be
+\begin{itemize}
+ \item $c_m = F\cdot cm^{-2}$
+ \item $V = V$
+ \item $a = cm$
+ \item $r_L = \Omega\cdot cm$
+\end{itemize}
+Then the units of each term in equation are $A\cdot cm^{-2}$.
+In practice, the units above are not used, for example distances are usually measured in $\mu m$ and areas in $cm^2$.
+But if we work term-by-term, scaling for these factors is manageable.
+
+A useful identity to use when performing the dimensional analysis relates capacitance and resistance
+\begin{equation*}
+ 1~F = 1~\Omega^{-1} \cdot s
+\end{equation*}
+
--
GitLab