An unconditionally stable numerical method for the Luo–Rudy 1 model used in simulations of defibrillation

Numerical simulations of defibrillation using the Bidomain model coupled to a model of membrane kinetics represent a serious numerical challenge. This is because very high voltages close to defibrillation electrodes demand that extreme time step restrictions be placed on standard numerical schemes, e.g. the forward Euler scheme. A common solution to this problem is to modify the cell model by simple if-tests applied to several equations and rate functions. These changes are motivated by numerical problems rather than physiology, and should therefore be avoided whenever possible.The purpose of this paper is to present a numerical scheme that handles the original model without modifications and which is unconditionally stable for the Luo–Rudy phase 1 model. This also shows that the cell model is mathematically well-behaved, even in the presence of very high voltages. Our theoretical results are illustrated by numerical computations.