\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
Eighth Mississippi State - UAB Conference on Differential Equations and
Computational Simulations.
{\em Electronic Journal of Differential Equations},
Conf. 19 (2010), pp. 189--196.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document} \setcounter{page}{189}
\title[\hfilneg EJDE-2010/Conf/19/\hfil Comparison of time stepping schemes]
{Comparison of time stepping schemes on the cable equation}
\author[C. Li, V. Alexiades \hfil EJDE/Conf/19 \hfilneg]
{Chuan Li, Vasilios Alexiades}
\address{Chuan Li \newline
Mathematics Department, University of Tennessee, Knoxville TN
37996, USA}
\email{li@math.utk.edu}
\address{Vasilios Alexiades \newline
Mathematics Department, University of Tennessee, Knoxville TN
37996, USA}
\email{alexiades@utk.edu}
\thanks{Published September 25, 2010.}
\subjclass[2000]{65M08, 35K57, 92C37}
\keywords{Explicit schemes; super time stepping; adaptive Runge Kutta;
\hfill\break\indent
Dufort Frankel; action potential; Luo-Rudy ionic models}
\begin{abstract}
Electrical propagation in excitable tissue, such as nerve
fibers and heart muscle, is described by a parabolic PDE
for the transmembrane voltage $V(x,t)$, known as the cable equation,
$$
\frac{1}{r_a}\frac{\partial^2V}{\partial x^2} =
C_m\frac{\partial V}{\partial t} + I_{\rm ion}(V,t)
+ I_{\rm stim}(t)
$$
where $r_a$ and $C_m$ are the axial resistance and membrane
capacitance. The source term $I_{\rm ion}$ represents the total ionic
current across the membrane, governed by the Hodgkin-Huxley or other
more complicated ionic models. $I_{\rm stim}(t)$ is an applied stimulus
current.
We compare the performance of various low and high order
time-stepping numerical schemes, including DuFort-Frankel
and adaptive Runge-Kutta, on the 1D cable equation.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\section{Biological Background}
It has been known since the 1940's that there is an electrical potential
difference between the inside and the outside of nerve cells.
However, researchers lacked methods to measure the membrane potential
directly at that time. In 1952, Hodgkin and Huxley published a series
of articles, based on experimental data obtained by others,
in which they established the first successful electrophysiological
model and unveiled the key properties of the ionic currents underlying
the nerve action potential \cite{HH:1952}. Inspired by their work and by
experimental developments, many other Hodgkin-Huxley type ionic models
have been created since then \cite{Foster:1993}.
Due to limited availability of human cardiomyocytes for experimental
research, most detailed electrophysiological models have been formulated
from animal cardiomyocytes. Luo-Rudy models \cite{LR1:1991, LR2a:1994,
LR2b:1994}, which were formulated for guinea pig ventricular cells,
are among the most widely used. In this paper, we focus on human
cardiac cells with the Luo-Rudy I (1991) ionic model \cite{cellML}.
\section{Mathematical Model}
A chain of cells, and a single cell with its equivalent circuit, are
represented in Figure \ref{fig:1} \cite{Shaw:1997} and
Figure \ref{fig:2} \cite{Jack:1975}, respectively.
The membrane is represented by discrete elements in parallel.
The membrane capacitance is ${C_m}$, the membrane resistance
is ${r_m}$ and the intracellular resistance
is represented by $r_a$.
\begin{figure}[htp]
\begin{center}
\includegraphics{fig1.jpg}
\end{center}
\caption{A typical single cell column model describing
discontinuous propagation}
\label{fig:1}
\end{figure}
\begin{figure}[htp]
\begin{center}
\includegraphics{fig2.jpg}
\end{center}
\caption{Cell with the axial and membrane currents and
its equivalent circuit}
\label{fig:2}
\end{figure}
The equation which describes the electrical behavior of this
system is a nonlinear parabolic equation
\cite{Keener:1998, Plonsey:2007},
known as the {\it cable equation}:
\begin{equation}
\frac{1}{r_a}\frac{\partial^2V}{\partial x^2} =
C_m\frac{\partial V}{\partial t} + I(V,t), \quad
\text{with } I(V,t) = I_{\rm ion}(V,t) + I_{\rm stim}(t),
\label{eqn1}
\end{equation}
where $V$ is the transmembrane voltage, $r_a$ and
$C_m$ are the axial resistance and membrane capacitance.
$I_{\rm ion}$ represents the total ionic current,
and $I_{\rm stim}(t)$ is the applied stimulus current.
In the Luo-Rudy 1991 model \cite{LR1:1991}, $I_{\rm ion}$
consists of several ionic currents generated by sodium,
potassium and calcium ions,
\begin{equation}
I_{\rm ion}(V,t) = I_{Na}(V) + I_{SI}(V) + I_{K1(T)}(V).
\label{eqn2}
\end{equation}
These three currents depend on seven activation and inactivation
``gates'': $m$, $h$, $j$, $d$, $f$, $X$, $Cai$,
which are governed by ODEs of the form
\begin{equation}
\frac{dg}{dt} = \alpha_{g} (V)(1-g) - \beta_{g}(V)g,
\quad g = m,h,j,d,f,X,Cai,
\label{eqn3}
\end{equation}
where the $\alpha$'s and $\beta$'s, taking values between 0 and 1,
are given by explicit formulas as functions of voltage $V$
\cite{LR1:1991}.
The stimulus current $I_{\rm stim}$, is the key to exciting the system.
In the heart, the stimulus is supplied by the Sino-Atrial Node.
Here we apply a single stimulus, of duration $1~\mu s$ and strength
$-200~\mu A/cm^2$, on a short segment $[0, 10~\mu m]$ at one end of
the cable.
The membrane capacitance $C_m$ is set at $1 ~\mu F/cm^2$ and the cable
radius $a$ is set at $4 ~\mu m$. The intracellular and extracellular
cytoplasmic resistivities are set to be $Ri(0) = 340~k\Omega cm$ and
$Ri(1)=1000~k\Omega cm$, respectively.
The corresponding resistances $r_a$ are computed by
\begin{equation}
r_a (k) = 2\pi a \frac{Ri(k)}{\pi a^2} = \frac{2Ri(k)}{a}
\quad for \quad k=0,1.
\label{eqn4}
\end{equation}
We assume the cable has insulted boundaries, and the system
starts from a steady state with initial values of
$ V_{\rm init} = -84.54816\, mV$, $m_{\rm init} = 0.00167$,
$h_{\rm init} = 0.9833$, $j_{\rm init} = 0.98952$,
$d_{\rm init} = 0.00298$, $f_{\rm init} = 0.99998$,
$X_{\rm init} = 0.0564$ and $Cai_{\rm init} = 0.0002$.
The model consists of the PDE \eqref{eqn1}, the seven
ODEs \eqref{eqn3} and the above initial setup.
See \cite{cellML} for more details.
\section{Numerical Schemes}
We discretize the cable into $M$ control volumes of uniform
length $\Delta x$ such that each cell contains multiple
control volumes.
Using standard Finite Volume discretization of the PDE \eqref{eqn1},
and applying equation \eqref{eqn3} on each control volume yields a
system of $8\,M$ ODEs,
\begin{equation}
\begin{gathered}
\frac{dV_k}{dt} =
\frac{1}{C_m} \Big[ \frac{F_{k-\frac{1}{2}}
- F_{k+\frac{1}{2}}}{\Delta x} - I(V_k, t_n) \Big] , \quad
I(V_k, t_n) = I_{\rm ion}(V_k, t_n) + I_{\rm stim}(t_n), \\
\frac{dg_k}{dt} = \alpha_{g_k}(V_k)(1-g_k) - \beta_{g_k}(V_k)g_k,
\quad k = 1,\dots ,M,
\end{gathered} \label{eqn8}
\end{equation}
where $V_k$ and $g_k = m_k$, $h_k$, $j_k$, $d_k$, $f_k$, $X_k$,
$Cai_k$ are the voltage and corresponding gates on the $k^{th}$
control volume, and $F_{k-\frac{1}{2}}$ and $F_{k+\frac{1}{2}}$
are the (diffusive) fluxes at the left and right faces, respectively,
\begin{equation}
F_{k-\frac{1}{2}} = \frac{1}{{r_a}_{k-\frac{1}{2}}} \frac{
V_{k-1} - V_{k}}{\Delta x}, \quad k = 2,\dots ,M
\label{eqn9}
\end{equation}
We apply the following numerical schemes to solve the
ODE system \eqref{eqn8}.
\subsection*{Super-Time-Stepping (STS) Scheme}
Super time-stepping is a simple way to accelerate explicit
schemes for parabolic problems \cite{Alexiades:1996}. One superstep
$\Delta T$ consists of $N$ substeps $\Delta \tau_1$,\dots ,$\Delta \tau_N$,
with optimal substeps $\Delta \tau_j$ given explicitly by
\begin{equation} %%%%%%%%%%%%%%% eq 7 %%%%%%%%%%%%%%%
\Delta \tau_{j} =
\Delta t_{\rm expl} \Big[ (-1+\nu)
\cos \Big(\frac{2j-1}{N}\frac{\pi}{2}\Big) + 1+\nu \Big]^{-1}
\quad j=1,\dots ,N,
\label{eqn5}
\end{equation}
where $\Delta t_{\rm expl}$ is the time step satisfying
the CFL stability
condition for the explicit scheme. Thus, we choose an integer $N$ and
a small damping parameter $\nu >0$, and instead of executing $N$
uniform steps $\Delta t_{\rm expl}$ we execute $N$ Chebyshev steps
$\Delta \tau_1,\dots ,\Delta \tau_N$. One can show the relation
\begin{equation}
\Delta T = \sum_{j=1}^{N} \Delta \tau_{j} =
\Delta t_\mathrm{expl}\frac{N}{2\sqrt{\nu}}
\Big[ \frac{(1+\sqrt{\nu })^{2N}-(1-\sqrt{\nu})^{2N}}
{(1+\sqrt{\nu })^{2N}+(1-\sqrt{\nu})^{2N}}\Big],
\label{eqn6}
\end{equation}
which yields
\begin{equation}
\Delta T \to N^2\Delta t_\mathrm{expl} \quad \text{as } \nu
\to 0.
\label{eqn7}
\end{equation}
Noting that $N$ explicit steps, each of length $\Delta t_{\rm expl}$,
cover time $N\Delta t_{\rm expl}$, we see that executing a superstep
consisting of $N$ substeps covers a time interval $N$ times longer
(when $\nu \approx 0$). Thus, superstepping is (up to) $N$ times
faster than the standard explicit scheme, at essentially the same
cost. This is where the speed-up comes from.
Note that the method ensures stability only at the end of each
superstep. Only the values at the end of each superstep approximate
the solution of the problem. STS reduces to plain Forward Euler by
setting parameters $N=1$, $\nu=0$.
In addition to speeding up the computation, the super-time-stepping
scheme is extremely simple to implement in any existing explicit code.
\subsection*{Adaptive Runge Kutta (RK) Schemes}
The package we chose to perform adaptive Runge Kutta methods is RKSUITE
from Netlib \cite{RKSuite:1991}. RKSUITE is a suite of codes based on
explicit Runge-Kutta methods for the numerical solution of the initial
value problem for a first order system of ordinary differential
equations. It supersedes some widely used codes.
Three adaptive methods, namely, RK23, RK45 and RK78 are provided in
this suite. Adaptive time step $\Delta t$ is controlled by two
parameters, relative tolerance $tol$ and threshold $thres$
provided by the user. To compare their performance, we set
$tol=10^{-3}$ and $thres=10^{-5}$ in all three methods for our tests.
\subsection*{DuFort-Frankel (DF) Scheme}
DuFort-Frankel is an explicit, 2-step, second order accurate in space
and time, theoretically unconditionally stable scheme
\cite{Morton:1994}. Applying centered finite difference in space and
first order Forward Euler in time on equation \eqref{eqn1} results in
\begin{equation}
V_k^{n+1} = V_k^n + \frac{\Delta t}{C_m \Delta x^2}
\Big( \frac{V_{k-1}^n - V_k^n}{{r_a}_{k-\frac{1}{2}}}
- \frac{V_{k}^n - V_{k+1}^n}{{r_a}_{k+\frac{1}{2}}} \Big)
+ \frac{\Delta t}{C_m} I(V_k^n, t_n).
\label{eqn10}
\end{equation}
Then $V_k^n$ is replaced by the average over two time steps
$(V_k^{n+1}+V_k^{n-1})/2$, producing a 2-step scheme.
To avoid small oscillations near the steady state, and keep the scheme
explicit, the average of voltage at previous two time steps is used to
evaluate the ionic current $I_{\rm ion}$,
\begin{equation}
\begin{aligned}
V_k^{n+1} &= \frac{1}{ 1+ \frac{\Delta t}{C_m \Delta x^2}
\big( \frac{1}{{r_a}_{k-\frac{1}{2}}} +
\frac{1}{{r_a}_{k+\frac{1}{2}}} \big) }
\Big[ \frac{2 \Delta t}{C_m \Delta
x^2}(\frac{1}{{r_a}_{k-\frac{1}{2}}}V_{k-1}^n +
\frac{1}{{r_a}_{k+\frac{1}{2}}}V_{k+1}^n) \\
&\quad
+ \Big(1 - \frac{\Delta t}{C_m \Delta x^2}
\big(\frac{1}{{r_a}_{k-\frac{1}{2}}} +
\frac{1}{{r_a}_{k+\frac{1}{2}}}\big)V_k^{n-1}\Big)
+ \frac{2 \Delta t}{C_m} I\big(\frac{V_k^n + V_k^{n-1}}{2}, t_n\big)
\Big].
\end{aligned} \label{eqn11}
\end{equation}
On the other hand, the ODEs \eqref{eqn3} for the gates are discretized
by forward Euler, and again evaluated at the average of the two
previous voltage values,
\begin{equation}
g_k^{n+1} = g_k^n + \Delta t \Big[
\alpha_{g_k^n} (\frac{V_k^n + V_k^{n-1}}{2})(1-g_k^n)
- \beta_{g_k^n}(\frac{V_k^n + V_k^{n-1}}{2})g_k^n \Big] .
\label{eqn12}
\end{equation}
Being aware of when the stimulus takes place, a time-step factor
$dtfac$ is introduced to speed up the computation here. That is,
$\Delta t_{\rm expl}$ is used in a small time period containing
the moment stimulus happens, and a larger time step
$\Delta t_{big}=dtfac*\Delta t_{\rm expl}$ is used elsewhere.
In the simulations, we used $dtfac=1$ and $dtfac=2$, denoting the
schemes as DF1 and DF2. They produce similar results.
\section{Numerical Simulations}
To compare the performance of the above numerical schemes, besides
execution time, the following biological quantities are significant:
\subsection*{Action Potential Duration (APD)} An action potential
is a transient alteration of the transmembrane potential
(voltage) generated by the activity of voltage-gated ion
channels embedded in the cell membrane.
The duration is determined by measuring how long the potential
$V$ at a location stays above a certain cut-off value.
In our computations, APD is
determined by setting a cut-off voltage 90\% of the initial
equilibrium voltage.
\subsection*{Propagation Speed (speed)} It measures how fast
the action potential propagates along the cable.
It is measured by the difference of the starting time of APDs
on the first and the last nodes.
\subsection*{Maximum voltage ($V_{\rm max}$) and maximum rate of
change ($dV/dt_{\rm max}$)}
These two quantities are calculated, for each scheme, at the nodes,
excluding those stimulated in the range $[0,{\rm stim}\_{\rm range}]$.
A C program has been written implementing the schemes. The numerical
experiments have been performed on a workstation equipped with dual AMD
Opteron 252 CPUs and 2GB memory, using the Portland Group pgCC compiler.
The resulting voltage history at
several equally separated nodes, and gates history at the first node
are shown in Figure \ref{fig:3}-\ref{fig:5}.
All schemes produce essentially identical plots.
Table 1 lists values obtained with each scheme on a cable of length
$5$ mm and one of length 10 mm. The last column shows the cost of each
computation (in seconds of CPU time).
\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig3.png}
% \resizebox*{0.65\textheight}{!}{{ \includegraphics{fig3.png}} }
\end{center}
\caption{Voltage history obtained by STS with $N=4$ and
$\nu = 0.07$ on a 5mm cable}
\label{fig:3}
\end{figure}
Based on our numerical experiments and results shown in Table 1,
we observe the following:
\begin{itemize}
\item All numerical schemes produce identical APD and very similar
propagation speeds.
\item The high order adaptive RK schemes are the most costly.
High cost of
evaluation of the source term $I_{\rm ion}$ is the key to this
phenomenon. To make RK family and other adaptive methods useful
in large size problems, we plan to test a ``library" method in
which $\alpha$'s and $\beta$'s in the ODEs are precomputed and
evaluated by interpolation as needed.
\item The first node, which receives full strength of the stimulus, has
significantly higher amplitude than the rest of the nodes which
receive stimulus via propagation. This indicates that
the ionic currents are highly sensitive to the change of
voltage. This high sensitivity restricts the choice of $dtfac$
used in the DF scheme. In fact, $dtfac=2$ is the best we can do.
Larger $dtfac$ makes the sodium gate ($m$) oscillate after
resting for a while and eventually blow up.
\item All schemes produce very similar $V_{\rm max}$ and $dV/dt_{\rm max}$ for
$5mm$ and $10mm$ cables and also for longer cables we tested.
\item On the basis of accuracy and efficiency (CPU time), STS and DF2
are the winners among the schemes tested.
\end{itemize}
\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig4.png}
% \resizebox*{0.65\textheight}{!}{ { \includegraphics{fig4.png}} }
\end{center}
\caption{Voltage history obtained by DF with $dtfac=2$
(DF2) on a 10mm cable}
\label{fig:4}
\end{figure}
\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig5.png}
% \includegraphics[width=0.96\textwidth,height=3.9in]{fig5.png}
\end{center}
\caption{Gates history at the 1st node obtained by RK45
on the 10mm cable}
\label{fig:5}
\end{figure}
\begin{table}
\caption{Comparison of timings on 1D cables with $\Delta x = 4 \mu m$
up to $t_{\rm max} = 2000ms$}
\label{table:1}
\begin{center}
\begin{tabular}{| l | c c c c r|}
\hline
\lower8pt\vbox to 22pt{}
scheme &APD[ms] &speed[cm/s] &$V_{\rm max}$[mV] &$dV/dt_{\rm max}$[V/s] &CPU[s] \\
\hline
&\multicolumn{5}{|l|}{\quad 5 mm cable} \\
FWD &371 &1.244 &40 &380 & ~515 \\
STS &371 &1.212 &42 &385 & ~284 \\
RK23 &371 &1.258 &41 &371 &2116 \\
RK45 &371 &1.258 &40 &370 &2815 \\
RK78 &371 &1.258 &41 &369 &5877 \\
DF1 &371 &1.246 &41 &379 & ~543 \\
DF2 &371 &1.234 &41 &379 & ~278 \\
\hline
&\multicolumn{5}{|l|}{\quad 10 mm cable} \\
FWD &371 &1.242 &40 &380 &1031 \\
STS &371 &1.210 &42 &385 & 573 \\
RK23 &371 &1.256 &41 &371 &4075 \\
RK45 &371 &1.256 &40 &370 &5630 \\
RK78 &371 &1.256 &41 &369 &10067 \\
DF1 &371 &1.244 &41 &379 &1075 \\
DF2 &371 &1.232 &41 &379 & 553 \\
\hline
\end{tabular}
\end{center}
\end{table}
\subsection*{Future Work}
We plan to implement and compare other numerical schemes and other
cardiac ionic models. Meanwhile, parallel codes for 1-, 2- and 3-
dimensional models are under development to speed up the computations.
Simulating cardiac arrhythmias is the goal.
\subsection*{Acknowledgments}
We thank Dr. Jack Buchanan of the University of Tennessee
Health Science Center in Memphis for providing biological parameters.
This work was supported by grant 1R21GM080698-01A1 from NIH.
\begin{thebibliography}{00}
\bibitem{Alexiades:1996} Vasilios Alexiades, Genevive Amiez,
and Pierre-Alain Gremaud,
\emph{Super-time-stepping acceleration of explicit schemes
for parabolic problems}, Commun. Num. Meth. Eng \textbf{12} (1996),
12--31.
\bibitem{RKSuite:1991}
RW~Brankin, I~Gladwell and LF~Shampine, \emph{RKSUITE},
\text{http://www.netlib.org/ode/rksuite/}.
\bibitem{cellML}
cellML, \emph{luo\_rudy\_1991\_version06},
\text{http://models.cellml.org/luo\_rudy\_1991\_version06}.
\bibitem{Foster:1993}
WR~Foster, LH~Unger and JS~Schwaber,
\emph{Significance of conductances in Hodgkin-Huxley models},
J. of Neuronphysiology \textbf{70} (1993).
\bibitem{HH:1952}
AL~Hodgkin and AF~Huxley,
\emph{A quantitative description of membrane current and its
application to conduction and excitation in nerve}, J.Physiol.
\textbf{117} (1952), 500--544.
\bibitem{Jack:1975}
JB~Jack, D~Noble and RW~Tsien, \emph{Electrical current flow
in excitable cells}, Oxford University Press, London, 1975.
\bibitem{Keener:1998}
James~Keener and James~Sneyd, \emph{Mathematical physiology},
Springer, 1998.
\bibitem{LR1:1991}
CH~Luo and Y~Rudy, \emph{A model of the ventricular cardiac action
potential: depolarization, repolarization, and their interaction},
Circ. Res. \textbf{68} (1991), 1501--1526.
\bibitem{LR2a:1994}
CH~Luo and Y~Rudy, \emph{A dynamic model of the cardiac ventricular
action potential. i. simulations of ionic currents and concentration
changes}, Circ. Res. \textbf{74} (1994), 1071–1096.
\bibitem{LR2b:1994}
CH~Luo and Y~Rudy, \emph{A dynamic model of the cardiac ventricular
action potential. ii. afterdepolarizations, triggered activity and
potentiation.}, Circ. Res. \textbf{74} (1994), 1097--1111.
\bibitem{Morton:1994}
DF~Mayers and KW~Morton, \emph{Numerical solution of partial
differential
equations}, 1st ed., Cambridge University Press, 1994.
\bibitem{Plonsey:2007}
Robert~Plonsey and Roger~C~Barr, \emph{Bioelectricity, a quantitative
approach}, 3rd ed., Springer, 2007.
\bibitem{Shaw:1997}
RM~Shaw and Y~Rudy, \emph{Electrophysiologic effects of acute myocardial
ischemia: a theoretical study of altered cell excitability and action
potential duration.}, Cardiovasc Res. \textbf{35} (1997), 181--183.
\end{thebibliography}
\end{document}