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