make correction to FVM formulation in docs

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}