在科學計算中，多維度心臟動作電位模擬非常耗時及需要佔用龐大記憶空間，因此運用高速率兼具準確的適應性演算法可以有效降低計算量。一般而言，使用高階Runge-Kutta (R-K)數值方法可獲得準確的結果，但造成計算時間冗長，所以開發準確兼備效率的適應性演算法是本研究的目標。心臟動作電位在去極化的相位中，其速度非常快速必需使用足夠小的時間步長以求得穏定計算；再者，位再極化相位中，則可使用適中的步長加速計算，然而在這些相位轉換區域間，當速度接近零時亦趨近極值，Rush-Larsen (R-L)及hybrid方法明顯選取過長的時間步長且無法辨識極值，因此時常造成計算中斷。針對這個缺點，本研究開發一種CCL(Chen-Chen-Luo)二次方適應性演算法在速度接近零而且趨近極值時，根據膜電位變化適應性地選擇小的時間步長，且在經過極植之後給予2倍增量逐步增加步長。此演算法，其由二次方程式組成配合計算動作電位的一次及二次微分進而適應性地擇取時間步長。為增進計算效率及穏定，應用時間步長限制器(time-step restriction)技巧將時間步長給予適時的抑制避免增量過度造成膜電位突然變化不穩定，以及設計另一極值定位器(extremum-locator)技術每當動作電位接近各峰值時，其可以發揮預測膜電壓極值發生之做用。CCL方法融合單細胞之羅-魯廸及歐-魯廸等心臟動作模型加速計算動作電位；並且成功植入Mono-domain模型應用Operator splitting演算法配合Crank-Nicolson及alternating direction implicit(ADI)方法建構1維及2維心臟細胞模擬平臺。數值結果顯示，CCL演算法所產生的模擬動作電位在最大峰值(Vmax)相較於高精度之RK4演算法約有0.06%的差異，當膜電位差(ΔV)被嚴格限制在0.01 mV時可獲得形廓流暢的模擬動作電位。針對0或1或2維度的動作電位模擬，CCL(0.1 mV)的誤差皆低於傳統的hybrid(0.05-0.2 mV)方法，其未顯著增加運算時間而且比RK4(0.01 ms)快約400倍。
本論文介紹二個臨床應用案例；首先，當CCL演算法應用在歐-魯廸心臟細胞模型，模擬Dofetilide藥物造成85% IKr抑制時，導致引發早發性後去極化(Early afterdepolarizations, EADs)動作電位之發生，況且其動作電位之速度接近零(極值)，CCL成功地辨識極值並選擇小步長且在經過峰值後可以逐漸增加時間步長變量加速計算，因此提昇計算的準確性及計算效率；其二，R-L方法應用在1維度的歐-魯廸心臟細胞模型時，在評估最佳時機緩解心搏過速的模擬實驗中，低電量應用在內皮層(Endocardial)及外皮層(Epicardial)心肌時，心搏過速皆可緩解；但是高電量應用在內皮層時，其會緩解心搏過速並引發自發性心律；反之，應用在外皮層則會造成心肌震顫(Ventricular Fibrillation) 。
總結而言，發表的CCL方法其計算效能是較準確於傳統的hybrid方法，而且計算時間也相當優異;未來引用Wenckebach perodicity 或Brugada sydromehqm 等生理現象表徵做為檢測驗證演算法穏定度之準則是成為持續精進的目標，此研究希望能解釋計算效能、準確性及穏定度之間的關係，研究人員也可藉此進一步了解數值方法是如何影響其動作電位模擬。在臨床電衝擊引發後效應模擬實驗中，其最佳效益緩解心搏過速是建議應用在心電圖之R波前實行。

In scientific computation, multi-dimensional cardiac action potential (AP) simulations require a considerably large amount of computation time and memory. To address this critical, a highly efficient and accurate adaptive numerical method that can reduce the amount of computation time is needed. In general, taking advantage of the accuracy of the Runge-Kutta (RK) method can obtain precise results, but this requires a lengthy computation time. Therefore, developing an accurate and efficient adaptive algorithm is paramount. During the depolarization phase of AP, its speed is ultra-fast and requires a very small time-step size to obtain a stable computation; furthermore, during the repolarization phase of AP it can use an adequate sized time-step to speed up the computation. However, during the transition zones between phases when the AP speed approaches zero, the Rush-Larsen (R-L) and hybrid methods cannot use small time-steps and identify the AP extremum value, which sometimes causes interruptions in computer simulations. In response to this drawback, the present study developed a quadratic formula for the CCL (Chen-Chen-Luo) adaptive method to increase the AP computation speed and maintain accuracy.
A quadratic CCL adaptive integration method is proposed to speed up the computation of the Luo-Rudy phase 1 (LR91) and O’Hara-Rudy ventricular AP models. Time-steps are adaptively chosen by solving a quadratic formula for the first and second derivative of the membrane action potential. To improve performance, the time-step restriction (tsr) technique is designed to limit the increase of time steps and thus avoid abrupt changes in membrane potential. Also, the extremum-locator (el) technique is used to predict the local extremum when approaching the peak AP amplitude. The CCL method was successfully coupled to the mono-domain model with cooperating Crank-Nicolson and Operator splitting methods to build up the computational platform for 1-Dimentional (1-D) cardiac fiber, and with cooperating alternating direction implicit (ADI) and Operator splitting methods to build up the computational platform for 2-Dimsensional (2-D) cardiac tissue simulations.
Numerical experiments show that the generated AP using the proposed method has a slight error of 0.06% at the peak (Vmax) compared to that of the fourth-order RK method (reference solution). The morphology of the AP is considerably improved when the membrane potential offset (ΔV) is strictly controlled by the small value of 0.01 mV. The CCL (0.1 mV) method is more accurate than the hybrid (0.05-0.2 mV) method, and is about 400% faster than that of the RK4 (0.01 ms) method in 0-D, 1-D and 2-D simulations. Overall, the proposed method is more accurate than the traditional hybrid method with an insignificant cost to computation time.
This paper introduces two clinical application cases. First, the CCL method is implemented for the O’Hara-Rudy model to simulate EADs induced by a drug (Dofetilide) effect due to the 85% IKr block. During EAD occurrences when the AP speed approaches zero, the CCL method chooses small time-step sizes, and, after passing the peak, the time steps are gradually increased two-fold to stabilize the computation. Following this, the R-L method is applied to a 1-D O’Hara-Rudy AP strand wedge model to evaluate the aftereffect of a Direct Current (DC) shock on endocardium and epicardium myocyte, in which a low-strength-DC shock is delivered to an endocardial or epicardial cell that can terminate the ventricular tachycardia (VT). By contrast, a high-strength-DC shock was delivered to an endocardial cell, enabling it to terminate VT and produce a spontaneous beat; however, a high-strength-DC shock delivered to an epicardial cell resulted in ventricular fibrillation (VF).
In conclusion, the computational performance of the proposed CCL numerical method is considerably more accurate than that of the traditional hybrid method. In the near future, the electrophysiology pattern of Wenckebach periodicity or Brugada syndrome (Brs) could be used to validate the computation stability for the CCL method and could potentially form a standard. This study may be of importance in explaining how the CCL method increases the computational performance, accuracy and stability, as well as in providing researchers with a better understanding of how numerical methods affect cardiac AP simulations. Additionally, according to the various DC shock-induced electroporation aftereffect simulations, the optimized timing of delivering a DC shock should be at the front of the R wave of the electrocardiogram.

[1] A. L. Hodgkin and A. F. Huxley, “A Quantitative Description of Membrane Current and Its Application to Conduction and Excitation in Nerve,” The Journal of physiology, vol. 117, pp. 500-544, 1952.
[2] A. L. Hodgkin and A. F. Huxley, “Movement of radioactive potassium and membrane current in a giant axon,” J Physiol, vol. 121, pp. 403-14, Aug 1953.
[3] J. W. Moore and F. Ramon, “Numerical Integration of Hodgkin and Huxley Equations for a Membrane Action Potential,” Journal of Theoretical Biology, vol. 45, pp. 249-273, 1974.
[4] S. Rush and H. Larsen, “Practical Algorithm for Solving Dynamic Membrane Equations,” IEEE Trans Biomed Eng, vol. 25, pp. 389-392, 1978.
[5] R. E. Mcallister, D. Noble, and R. W. Tsien, “Reconstruction of Electrical-Activity of Cardiac Purkinje-Fibers,” Journal of Physiology-London, vol. 251, pp. 1-59, 1975.
[6] B. Victorri, A. Vinet, F. A. Roberge, and J. P. Drouhard, “Numerical-Integration in the Reconstruction of Cardiac Action-Potentials Using Hodgkin-Huxley-Type Models,” Computers and Biomedical Research, vol. 18, pp. 10-23, 1985.
[7] J. Sundnes, R. Artebrant, O. Skavhaug, and A. Tveito, “A Second-Order Algorithm for Solving Dynamic Cell Membrane Equations,” IEEE Trans Biomed Eng, vol. 56, pp. 2546-2548, Oct 2009.
[8] S. Linge, J. Sundnes, M. Hanslien, G. T. Lines, and A. Tveito, “Numerical solution of the bidomain equations,” Philosophical Transactions of the Royal Society a-Mathematical Physical and Engineering Sciences, vol. 367, pp. 1931-1950, May 28 2009.
[9] W. J. Ying, D. J. Rose, and C. S. Henriquez, “Efficient Fully Implicit Time Integration Methods for Modeling Cardiac Dynamics,” IEEE Trans Biomed Eng, vol. 55, pp. 2701-2711, 2008.
[10] R. J. Spiteri and R. C. Dean, “On the performance of an implicit-explicit Runge-Kutta method in models of cardiac electrical activity,” IEEE Trans Biomed Eng, vol. 55, pp. 1488-1495, May 2008.
[11] J. P. Whiteley, “An efficient numerical technique for the solution of' the monodomain and bidomain equations,” IEEE Trans Biomed Eng, vol. 53, pp. 2139-2147, Nov 2006.
[12] R. J. Spiteri and M. C. MacLachlan, “An efficient non-standard finite difference scheme for an ionic model of cardiac action potentials,” Journal of Difference Equations and Applications, vol. 9, pp. 1069-1081, Sep 2003.
[13] Z. L. Qu and A. Garfinkel, “An advanced algorithm for solving partial differential equation in cardiac conduction,” IEEE Trans Biomed Eng, vol. 46, pp. 1166-1168, Sep 1999.
[14] J. P. Keener and K. Bogar, “A numerical method for the solution of the bidomain equations in cardiac tissue,” Chaos, vol. 8, pp. 234-241, Mar 1998.
[15] D. Noble, “A modification of the Hodgkin--Huxley equations applicable to Purkinje fibre action and pace-maker potentials,” The Journal of physiology, vol. 160, pp. 317-52, Feb 1962.
[16] C. H. Luo and Y. Rudy, “A Model of the Ventricular Cardiac Action-Potential - Depolarization, Repolarization, and Their Interaction,” Circulation Research, vol. 68, pp. 1501-1526, Jun 1991.
[17] S. C. C. a. R. P. Canale, “Numerical Methods for Engineers With Software and Programmimg Applications,” New York, USA: McGraw-Hill, 2002.
[18] K. Gima and Y. Rudy, “Ionic Current Basis of Electrocardiographic Waveforms,” Circulation Research, vol. 90, pp. 889-896, May 3, 2002 2002.
[19] T. O'Hara, L. Virág, A. Varró, and Y. Rudy, “Simulation of the Undiseased Human Cardiac Ventricular Action Potential: Model Formulation and Experimental Validation,” PLoS Comput Biol, vol. 7, p. e1002061, 2011.
[20] M. Josephson and H. J. Wellens, “Implantable defibrillators and sudden cardiac death,” Circulation, vol. 109, pp. 2685-91, Jun 8 2004.
[21] B. F. Uretsky and R. G. Sheahan, “Primary prevention of sudden cardiac death in heart failure: will the solution be shocking?,” J Am Coll Cardiol, vol. 30, pp. 1589-97, Dec 1997.
[22] W. G. Stevenson, L. W. Stevenson, H. R. Middlekauff, and L. A. Saxon, “Sudden death prevention in patients with advanced ventricular dysfunction,” Circulation, vol. 88, pp. 2953-61, Dec 1993.
[23] A. Bayes de Luna, P. Coumel, and J. F. Leclercq, “Ambulatory sudden cardiac death: mechanisms of production of fatal arrhythmia on the basis of data from 157 cases,” Am Heart J, vol. 117, pp. 151-9, Jan 1989.
[24] C. S. Fox, J. C. Evans, M. G. Larson, W. B. Kannel, and D. Levy, “Temporal trends in coronary heart disease mortality and sudden cardiac death from 1950 to 1999: the Framingham Heart Study,” Circulation, vol. 110, pp. 522-7, Aug 3 2004.
[25] T. J. Bunch, R. D. White, P. A. Friedman, T. E. Kottke, L. A. Wu, and D. L. Packer, “Trends in treated ventricular fibrillation out-of-hospital cardiac arrest: a 17-year population-based study,” Heart Rhythm, vol. 1, pp. 255-9, Sep 2004.
[26] A. Kadish, A. Dyer, J. P. Daubert, R. Quigg, N. A. Estes, K. P. Anderson, H. Calkins, D. Hoch, J. Goldberger, A. Shalaby, W. E. Sanders, A. Schaechter, and J. H. Levine, “Prophylactic defibrillator implantation in patients with nonischemic dilated cardiomyopathy,” N Engl J Med, vol. 350, pp. 2151-8, May 20 2004.
[27] A. J. Moss, W. Zareba, W. J. Hall, H. Klein, D. J. Wilber, D. S. Cannom, J. P. Daubert, S. L. Higgins, M. W. Brown, and M. L. Andrews, “Prophylactic implantation of a defibrillator in patients with myocardial infarction and reduced ejection fraction,” N Engl J Med, vol. 346, pp. 877-83, Mar 21 2002.
[28] S. B. Eysmann, F. E. Marchlinski, A. E. Buxton, and M. E. Josephson, “Electrocardiographic changes after cardioversion of ventricular arrhythmias,” Circulation, vol. 73, pp. 73-81, Jan 1986.
[29] J. N. Patton, J. D. Allen, and J. F. Pantridge, “The effects of shock energy, propranolol, and verapamil on cardiac damage caused by transthoracic countershock,” Circulation, vol. 69, pp. 357-68, Feb 1984.
[30] B. B. Lerman, J. L. Weiss, B. H. Bulkley, L. C. Becker, and M. L. Weisfeldt, “Myocardial injury and induction of arrhythmia by direct current shock delivered via endocardial catheters in dogs,” Circulation, vol. 69, pp. 1006-12, May 1984.
[31] J. L. Jones and R. E. Jones, “Postshock arrhythmias--a possible cause of unsuccessful defibrillation,” Crit Care Med, vol. 8, pp. 167-71, Mar 1980.
[32] T. O'Hara, L. Virag, A. Varro, and Y. Rudy, “Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation,” PLoS Comput Biol, vol. 7, p. e1002061, May 2011.
[33] K. Ohuchi, Y. Fukui, I. Sakuma, N. Shibata, H. Honjo, and I. Kodama, “A dynamic action potential model analysis of shock-induced aftereffects in ventricular muscle by reversible breakdown of cell membrane,” IEEE Trans Biomed Eng, vol. 49, pp. 18-30, Jan 2002.
[34] I. Sakuma, T. Haraguchi, K. Ohuchi, Y. Fukui, I. Kodama, J. Toyama, N. Shibata, and S. Hosoda, “A model analysis of aftereffects of high-intensity DC stimulation on action potential of ventricular muscle,” IEEE Trans Biomed Eng, vol. 45, pp. 258-267, 1998.
[35] D. Guo, Q. Liu, T. Liu, G. Elliott, M. Gingras, P. R. Kowey, and G. X. Yan, “Electrophysiological properties of HBI-3000: a new antiarrhythmic agent with multiple-channel blocking properties in human ventricular myocytes,” J Cardiovasc Pharmacol, vol. 57, pp. 79-85, 2011.
[36] C. Luo and Y. Rudy, “A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes,” Circulation Research, vol. 74, pp. 1071-1096, June 1, 1994.
[37] A. Al-Khadra, V. Nikolski, and I. R. Efimov, “The Role of Electroporation in Defibrillation,” Circulation Research, vol. 87, pp. 797-804, October 27, 2000.
[38] S. C. Kim, A. Vasanji, I. R. Efimov, and Y. Cheng, “Spatial Distribution and Extent of Electroporation by Strong Internal Shock in Intact Structurally Normal and Chronically Infarcted Rabbit Hearts,” Journal of Cardiovascular Electrophysiology, vol. 19, pp. 1080-1089, 2008.
[39] T. Ashihara and N. A. Trayanova, “Asymmetry in membrane responses to electric shocks: insights from bidomain simulations,” Biophys J, vol. 87, pp. 2271-82, Oct 2004.
[40] C. M. Ripplinger, V. I. Krinsky, V. P. Nikolski, and I. R. Efimov, “Mechanisms of unpinning and termination of ventricular tachycardia,” American Journal of Physiology - Heart and Circulatory Physiology, vol. 291, pp. H184-H192, July 2006.
[41] J. Zeng and Y. Rudy, “Early afterdepolarizations in cardiac myocytes: mechanism and rate dependence,” Biophysical Journal, vol. 68, pp. 949-964, 1995.
[42] W. Yuan and D. M. Bers, “Ca-dependent facilitation of cardiac Ca current is due to Ca-calmodulin-dependent protein kinase,” American Journal of Physiology - Heart and Circulatory Physiology, vol. 267, pp. H982-H993, September 1, 1994.
[43] E. Soltysinska, S.-P. Olesen, T. Christ, E. Wettwer, A. Varró, M. Grunnet, and T. Jespersen, “Transmural expression of ion channels and transporters in human nondiseased and end-stage failing hearts,” Pflügers Archiv European Journal of Physiology, vol. 459, pp. 11-23, 2009.
[44] N. Gaborit, S. Le Bouter, V. Szuts, A. Varro, D. Escande, S. Nattel, and S. Demolombe, “Regional and tissue specific transcript signatures of ion channel genes in the non-diseased human heart,” The Journal of Physiology, vol. 582, pp. 675-693, July 15, 2007.
[45] G. Szabó, N. Szentandrássy, T. Bíró, B. I. Tóth, G. Czifra, J. Magyar, T. Bányász, A. Varró, L. Kovács, and P. P. Nánási, “Asymmetrical distribution of ion channels in canine and human left-ventricular wall: epicardium versus midmyocardium,” Pflügers Archiv European Journal of Physiology, vol. 450, pp. 307-316, 2005.
[46] M. Courtemanche, “Complex spiral wave dynamics in a spatially distributed ionic model of cardiac electrical activity,” Chaos, vol. 6, pp. 579-600, 1996.
[47] D. Noble, A. Garny, and P. J. Noble, “How the Hodgkin-Huxley equations inspired the Cardiac Physiome Project,” J Physiol, vol. 590, pp. 2613-28, 2012.