From b0a010db10cabe3724d4d131136bdecfa762355b Mon Sep 17 00:00:00 2001
From: bcumming <bcumming@cscs.ch>
Date: Wed, 20 Jan 2016 10:33:45 +0100
Subject: [PATCH] make correction to FVM formulation in docs
---
docs/formulation.tex | 111 ++++++++++++++++++++-----------------------
1 file changed, 52 insertions(+), 59 deletions(-)
diff --git a/docs/formulation.tex b/docs/formulation.tex
index bf8f2fe4..3275b658 100644
--- a/docs/formulation.tex
+++ b/docs/formulation.tex
@@ -6,6 +6,7 @@ See \cite{lindsay_2004} for a detailed derivation of the cable equation, and ext
The one-dimensional cable equation introduced later in equations~\eq{eq:cable} and~\eq{eq:cable_balance} is based on the following expression in three dimensions (based on Maxwell's equations adapted for neurological modelling)
\begin{equation}
\nabla \cdot \vv{J} = 0,
+ \label{eq:J}
\end{equation}
where $\vv{J}$ is current density (units $A/m^2$).
Current density is in turn defined in terms of electric field $\vv{E}$ (units $V/m$)
@@ -53,27 +54,38 @@ where
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.
-The PDE in (\ref{eq:cable}) is derived from the following mass balance expression for a cable segment $i$:
-%\begin{align}
- %\int_{\Omega_i}{c_m \pder{V}{t} } \deriv{v} =
- %& + \sum_{j\in\mathcal{N}_i} {\int_{\Gamma_{i,j}} \frac{1}{r_L}\pder{V}{x} n_{i,j} \deriv{s} } \nonumber \\
- %& + \int_{\Gamma_{ext}} {(i_m - i_e)} \deriv{s}
- %\label{eq:cable_balance}
-%\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_{i}} {q_i} \deriv{s}
- \label{eq:cable_balance}
+The PDE in (\ref{eq:cable}) is derived by integrating~\eq{eq:J} over the volume of segment $i$:
+\begin{equation*}
+ \int_{\Omega_i}{\nabla \cdot \vv{J} } \deriv{v} = 0,
+\end{equation*}
+Then applying the divergence theorem to turn the volume integral on the lhs into a surface integral
+\begin{equation}
+ \sum_{j\in\mathcal{N}_i} {\int_{\Gamma_{i,j}} J_{i,j} \deriv{s} }
+ + \int_{\Gamma_{i}} {J_m} \deriv{s} = 0
+ \label{eq:cable_balance_intermediate}
\end{equation}
where
\begin{itemize}
\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 (where $q_i>0$ implies flux out of the cell).
+ \item $J_{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 the set $\mathcal{N}_i$ is the set of segments that are neighbours of $\Omega_i$
\end{itemize}
+The transmembrane current density $J_m=\vv{J}\cdot\vv{n}$ is a function of the membrane potential
+\begin{equation}
+ J_m = c_m\pder{V}{t} + i_m - i_e,
+ \label{eq:Jm}
+\end{equation}
+which has contributions from the ion channels and synapses ($i_m$), electrodes ($i_e$) and capacitive current due to polarization of the membrane whose bi-layer lipid structure causes it to behave locally like a parallel plate capacitor.
+
+Substituting~\eq{eq:Jm} into~\eq{eq:cable_balance_intermediate} and rearanging gives
+\begin{equation}
+ \int_{\Gamma_{i}} {c_m\pder{V}{t}} \deriv{s}
+ = - \sum_{j\in\mathcal{N}_i} {\int_{\Gamma_{i,j}} J_{i,j} \deriv{s} }
+ - \int_{\Gamma_{i}} {(i_m - i_e)} \deriv{s}
+ \label{eq:cable_balance}
+\end{equation}
+
The surface of the cable segment is sub-divided into the internal and external surfaces.
The external surface $\Gamma_{i}$ is the cell membrane at the interface between the extra-cellular and intra-cellular regions.
@@ -105,50 +117,42 @@ The finite volume method is a natural choice for the solution of the conservatio
%-------------------------------------------------------------------------------
\subsubsection{Temporal derivative}
%-------------------------------------------------------------------------------
-We proceed by defining the \emph{volume average} of a quantity $\varphi$ as follows:
-\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 $\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$
+The integral on the lhs of~\eq{eq:cable_balance} can be approximated by assuming that the average transmembrane potential $V$ in $\Omega_i$ is equal to the potential $V_i$ defined at the centre of the segment:
\begin{equation}
- \int_{\Omega}{c_m \pder{V}{t} } \deriv{v} = \Delta_i c_m \pder{\bar{V}_i}{t}.
+ \int_{\Gamma_i}{c_m \pder{V}{t} } \deriv{v} \approx \sigma_i c_m \pder{V_i}{t},
\label{eq:dvdt}
\end{equation}
-
-In the FV formulation the voltage $V_i$ at the node $x_i$ is equal to the volume average, i.e. $V_i=\bar{V}_i$.
-This is effectively treats voltage as a piecewise continuous funtion, with discontinuities at the boundary between adjacent segements.
+where $\sigma_i$ is the surface area of the membrane potential.
%-------------------------------------------------------------------------------
\subsubsection{Intra-cellular flux}
%-------------------------------------------------------------------------------
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}} { q_{i,j} \deriv{s} } }.
+ \sum_{j\in\mathcal{N}_i} { \int_{\Gamma_{i,j}} { J_{i,j} \deriv{s} } }.
\end{equation}
where the flux per unit area from compartment $i$ to compartment $j$ is
\begin{align}
- q_{i,j} = - \frac{1}{r_L}\pder{V}{x} n_{i,j}.
- \label{eq:q_ij}
+ J_{i,j} = - \frac{1}{r_L}\pder{V}{x} n_{i,j}.
+ \label{eq:J_ij_exact}
\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
+Using this approximation for the derivative, the flux over the surface in~\eq{eq:J_ij_exact} 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}
+ J_{i,j} \approx \frac{1}{r_L}\frac{V_i - V_j}{\Delta x_{i,j}}.
+ \label{eq:J_ij_intermediate}
\end{align}
-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
+The terms inside the integral in equation~\eq{eq:J_ij_intermediate} are constant everywhere on the surface $\Gamma_{i,j}$, so the integral becomes
\begin{align}
- q_{i,j} &= \int_{\Gamma_{i,j}} \frac{1}{r_L}\frac{V_i-V_j}{\Delta x_{i,j}} \deriv{s} \nonumber \\
- &= \frac{1}{r_L}\frac{V_i-V_j}{\Delta x_{i,j}} \int_{\Gamma_{i,j}} 1 \deriv{s} \nonumber \\
- &= \frac{1}{r_L}\frac{V_i-V_j}{\Delta x_{i,j}} \sigma_{i,j} \nonumber \\
- &= \frac{\pi a_{i,j}^2}{r_L \Delta x_{i,j}} (V_i-V_j)
- \label{eq:q_ij}
+ J_{i,j} &\approx \int_{\Gamma_{i,j}} \frac{1}{r_L}\frac{V_i-V_j}{\Delta x_{i,j}} \deriv{s} \nonumber \\
+ &= \frac{1}{r_L \Delta x_{i,j}}(V_i-V_j) \int_{\Gamma_{i,j}} 1 \deriv{s} \nonumber \\
+ &= \frac{\sigma_{ij}}{r_L \Delta x_{i,j}}(V_i-V_j) \nonumber \\
+ \label{eq:J_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_{i,j}$.
@@ -157,7 +161,7 @@ Some symmetries
\item $\sigma_{i,j}=\sigma_{j,i}$ : surface area of $\Gamma_{i,j}$
\item $\Delta x_{i,j}=\Delta x_{j,i}$ : distance between $x_i$ and $x_j$
\item $n_{i,j}=-n_{j,i}$ : surface ``norm''/orientation
- \item $q_{i,j}=n_{j,i}q_{i,j}=-q_{j,i}$ : charge flux over $\Gamma_{i,j}$
+ \item $J_{i,j}=n_{j,i}\cdot J_{i,j}=-J_{j,i}$ : charge flux over $\Gamma_{i,j}$
\end{itemize}
%-------------------------------------------------------------------------------
@@ -165,17 +169,17 @@ Some symmetries
%-------------------------------------------------------------------------------
The final term in~\eq{eq:cable_balance} with an integral is the cell membrane flux contribution
\begin{equation}
- q_{i}^{\text{m}} = \int_{\Gamma_{ext}} {(i_m - i_e)} \deriv{s},
+ \int_{\Gamma_{ext}} {(i_m - i_e)} \deriv{s},
\end{equation}
where the current $i_m$ is due to ion channel and synapses, and $i_e$ is any artificial electrode current.
The $i_m$ term is dependent on the potential difference over the cell membrane $V_i$.
The current terms are an average per unit area, therefore the total flux
-\begin{align}
- q_{i}^{\text{m}}
- & = (i_m(V_i) - i_e(x_i))\int_{\Gamma_{ext}} {1} \deriv{s} \nonumber \\
- & = \sigma_i\cdot(i_m(V_i) - i_e(x_i)),
- \label{eq:q_im}
-\end{align}
+\begin{equation}
+ \int_{\Gamma_{ext}} {(i_m - i_e)} \deriv{s}
+ \approx
+ \sigma_i(i_m(V_i) - i_e(x_i)),
+ \label{eq:J_im}
+\end{equation}
where $\sigma_i$ is the surface area the of the exterior of the cable segment, i.e. the surface corresponding to the cell membrane.
Each cable segment is a conical frustrum, as illustrated in \fig{fig:segment}.
@@ -190,18 +194,12 @@ 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}
- \Delta_i c_m \dder{V_i}{t}
+ \sigma_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
-%\begin{itemize}
-% \item $\sigma_{i,j}=\pi a_{i,j}^2$ is the area of the surface between two adjacent segments $i$ and $j$.
-% \item $\sigma_{i}=\pi(a_{i,\ell} + a_{i,r})+\sqrt{\Delta x_i^2 + (a_{i,\ell} - a_{i,r})^2}.$ is the lateral area of the conical frustrum describing segment $i$.
-% \item $\Delta_{i}=\frac{\pi\Delta x_i}{2} \left( a_{i,l}^2 + a_{i,r}^2 \right)$ is the volume of the segment $\Omega_i$.
-%\end{itemize}
\begin{equation}
\sigma_{i,j} = \pi a_{i,j}^2
\label{eq:sigma_ij}
@@ -211,23 +209,18 @@ is the area of the surface between two adjacent segments $i$ and $j$, and
\sigma_{i} = \pi(a_{i,\ell} + a_{i,r}) \sqrt{\Delta x_i^2 + (a_{i,\ell} - a_{i,r})^2},
\label{eq:sigma_i}
\end{equation}
-is the lateral area of the conical frustrum describing segment $i$, and
-\begin{equation}
- \Delta_{i} = \frac{\pi\Delta x_i}{2} \left( a_{i,l}^2 + a_{i,r}^2 \right)
- \label{eq:delta_i}
-\end{equation}
-is the volume of the segment $\Omega_i$.
+is the lateral area of the conical frustrum describing segment $i$.
%-------------------------------------------------------------------------------
\subsubsection{Handling branches}
%-------------------------------------------------------------------------------
-The value of the lateral area and volume, $\sigma_i$ and $\Delta_i$ in~\eq{eq:sigma_i} and~\eq{eq:delta_i} respetively, must include contributions from each branch at branch points.
+The value of the lateral area $\sigma_i$ in~\eq{eq:sigma_i} is the sum of the areaof each branch at branch points.
\todo{a picture of a branching point to illustrate}
\todo{a picture of a soma to illustrate the ball and stick model with a sphere for the soma and sticks for the dendrites branching off the soma.}
\begin{equation}
- \sigma_i = \sum_{j\in\mathcal{N}_i} {}
+ \sigma_i = \sum_{j\in\mathcal{N}_i} {\dots}
\end{equation}
--
GitLab