Exponential integrators for a Markov chain model of the fast sodium channel of cardiomyocytes

TL;DR

This paper introduces two exponential integrator methods, Matrix Rush-Larsen and hybrid operator splitting, to improve numerical stability and efficiency in simulating stiff Markov chain models of cardiac ion channels.

Contribution

The paper develops and tests two novel exponential integrator methods that enable larger time steps in Markov chain models of ion channels, reducing computational cost while maintaining accuracy.

Findings

Both methods allow longer time steps without instability.

They are comparable to forward Euler in accuracy and cost.

Matrix Rush-Larsen is more generalizable.

Abstract

The modern Markov chain models of ionic channels in excitable membranes are numerically stiff. The popular numerical methods for these models require very small time steps to ensure stability. Our objective is to formulate and test two methods addressing this issue, so that the timestep can be chosen based on accuracy rather than stability. Both proposed methods extend Rush-Larsen technique, which was originally developed to Hogdkin-Huxley type gate models. One method, "Matrix Rush-Larsen" (MRL) uses a matrix reformulation of the Rush-Larsen scheme, where the matrix exponentials are calculated using precomputed tables of eigenvalues and eigenvectors. The other, "hybrid operator splitting" (HOS) method exploits asymptotic properties of a particular Markov chain model, allowing explicit analytical expressions for the substeps. We test both methods on the Clancy and Rudy (2002) INa Markov chain model. With precomputed tables for functions of the transmembrane voltage, both methods are comparable to the forward Euler method in accuracy and computational cost, but allow longer time steps without numerical instability. We conclude that both methods are of practical interest. MRL requires more computations than HOS, but is formulated in general terms which can be readily extended to other Markov Chain channel models, whereas the utility of HOS depends on the asymptotic properties of a particular model. The significance of the methods is that they allow a considerable speed-up of large-scale computations of cardiac excitation models by increasing the time step, while maintaining acceptable accuracy and preserving numerical stability.