Exponential integrators for a Markov chain model
of the fast sodium channel of cardiomyocytes
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) 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.
Mathematical models are an essential part of the modern cardiac electrophysiology. They are used for hypothesis testing in research and as a guide for clinical decision. A typical definition of such a model is a high-dimensional (tens of equations) system of ordinary differential equations per excitable unit. Detailed simulations of the heart involve solving such systems for each of millions of cells placed in a mesh representing the cardiac tissue. Such large-scale models can be computationally extremely expensive, hence significant efforts are directed to develop efficient numerical methods for solving such systems.
A typical cardiac excitation model is centered around the Kirchhoff circuit law which gives
where is the cell membrane’s capacitance, is the transmembrane potential difference, and , , are currents through ion-specific channels. The currents, in turn, are determined by the Ohm’s law,
where is the ion-specific electromotive force, depending on the ionic concentrations via the Nerst equation, is the total conductance of channels of type when they are all open, and is the probability of those channels to be open.
The components of the vector are intra- and extra-cellular ionic concentrations, which change in time in the obvious way in accordance with the ionic fluxes and the corresponding volumes; some concentrations in some models are assumed constant. The dynamics of the open probabilities is much more nontrivial, as it reflects the conformation dynamics of the proteins constituting the ion channels .
The classical description of these dynamics, going back to Hodgkin and Huxley (1952) , has the form
with a popular, although different from the original Hodgkin and Huxley’s, interpretation that the set corresponds the subunits of the channel, called “gates”. These subunits are assumed statistically independent, each of them can be either in an “open” or a “closed” state, and the channel is open if and only if each of the subunits is open. Variables then are open probabilities of the gates, and their dynamics are described by
where are opening rates and are closing rates.
For instance, the original Hodgkin-Huxley description of the fast sodium current () channel uses gates, three of which, called m-gates, have identical opening and closing rates, and the fourth, called h-gate, has rates and , hence for this case we have
A more recent approach is modelling the probabilities of the channel molecules, as a whole, to be in specific conformation states, without the restricting assumptions of statistically independent subunits and only two states for any subunit. This gives generic continuous time Markov chain (MC) models. Let the probability of the ’th channel to be in the ’th state at time be (“state occupancy”), , , and all such probabilities be considered components of the state (column-) vector . Let be the relative permeability of the state , then we have
Typically, where is the “open state”. The time evolution of the state vector is described by the system of linear ODEs, known in particular as Kolmogorov (forward) equations, or master equation, of the form
where the non-diagonal components of the matrix are the transition rates (TR) between the states, and the diagonal components are defined by the condition and consequently sum of any column of should be zero.
The ODE system for cellular membrane can be solved on a computer using standard numerical solvers. A typical solver iteratively computes the states of the system using time-stepping algorithms, that is computing the states at times . The size of the time step is inversely proportional to the computational cost, measured as CPU time required for the computation. Increasing the time step is a straightforward way of reducing the computational cost.
The maximal acceptable time step is limited by considerations of accuracy and stability (see e.g. [14, Sections 5.10, 5.11]). Whereas the former is “relative” in that it depends on the aims of the research, the latter has a more “absolute” character in that if stability conditions are not satisfied, the solution is unusable for any purpose. Typically, when the time step exceeds the stability limit, the numerical solution is characterized by wild oscillations around the exact solution, and quite often will lead to numerical overflow.
Simple explicit solvers suffer from instabilities the most, and implicit, stable methods, applicable to generic systems of ODEs, are complicated and often costly. The motivation for our research was that taking into account the specific properties of the problem can offer some advantages. Specifically, we have in mind two distinct considerations.
One consideration is that the TRs can range through several orders of magnitude, and some of them can be much faster than other processes described by the excitable cell model. This split of the speeds of the variables suggest a possibility to exploit asymptotic methods.
The other consideration is the linearity of the system (9). Here we are inspired by the example of the exponential integrator algorithm developed by Rush and Larsen in 1978 . It is based on the assumption that the transmembrane voltage, on which the TRs in the gate model (4) depend, changes only slightly during one time step . So during one time step, the TRs can be approximated by constants, and the equation (4) is then solved analytically. The solution can be conveniently defined in terms of the “steady state” and the “time constant” at a given potential presumed constant:
The Rush-Larsen (RL) scheme is easy to implement, gives good results and is very popular in computational cardiac electrophysiology. Its stability and approximation properties have been extensively discussed in literature, including its relation to general exponential integrators family, its extension beyond gating variables by linearization, and improving its approximation properties, see e.g. [16, 17, 18]. However it is designed for a single ODE and is not immediately applicable for MC models which are systems of coupled ODEs. And yet MC models are known to suffer from severe numerical instability issues, just as, or even more than, the gate models (Fig. 1). The classical techniques for numerical solution of continuous-time MC models involve finding the eigenvalues and associated eigenvectors of the transition matrix. Direct implementation of this approach to very large MCs is problematic, see e.g. . However the MCs describing ionic channels are relatively small so the direct approach is feasible.
In this paper, we discuss two methods for numerical solution of MC models based on these two considerations.
To test the suggested numerical methods we have chosen the MC model of the channel by Clancy and Rudy  (Fig. 2), which is one of the most popular MC models. We used the formulation of the MC model and the whole cell model into which it was incorporated, as implemented in the authors’ code kindly provided by C.E. Clancy. It most closely corresponds to the Luo-Rudy model  with modifications described in [22, 23], and some further minor differences. For the sake of reproducibility of our results, we describe the whole model in the supplementary material, highlighting all the differences from the published models that we have detected. For the same purpose, we put a simplified version of the C code we used in the simulations described below in the supplementary materials.
For convenience, we changed the notation for the MC states and TRs. The states were named in alphabetical order, starting with for the open state, in a clockwise direction as in Fig. 2. See Tab. I and Tab. II for the correspondence with the original notation. The model contains 9 interconnected states. The state represents the conformation of the ion channel that allows the flow of ions between the intracellular and extracellular environment. The remaining states (, , , , , , and ) represent non-conductive conformations of the channel, so we can say that for this model , where . There are 11 possible bidirectional transitions between states, but some of the corresponding 22 TRs are described by identical functions, so there are only 14 distinct TR definitions. We denote the TRs by with a subscript showing the direction of the transition, e.g. is the transition rate from state into state . See Tab. II for the link with the original notations.
The TRs are shown on Fig. 3 as functions of the transmembrane potential in a physiologically relevant range. The values of TRs vary across several orders of magnitude, from to . Some of the TRs are high at the lower potentials, some are fast at higher potentials, and some are uniformly low.
The conductive (open) state is the only state that has immediate effect on the current. The remaining 8 states of the model can affect the current only indirectly by transitions to the open state . The time evolution of a generic state occupancy state is described by a differential equation of the form
where is the set of all the states interconnected with state , which can be readily found from the diagram. For example, the occupancy of the open state, , is described by the following equation:
By taking the sum of the equations (11) for all , one can see that the sum of the right-hand sides equals to zero, and therefore the system observes a states conservation law, which is consistent with the definition of as probabilities, implying . This is of course a generic property of a continuous Markov chain. So the differential equations in the model are not independent, which creates a possibility of reducing the number of equations from 9 to 8, by computing one of the occupancies through the conservation law rather than from its differential equation. However the computational gain from this is insignificant, and instead we used any deviations from the conservation law as an indicator of the accuracy of the computations.
Ii-B Numerical Methods
Ii-B1 Forward Euler
The standard forward Euler (FE) method is the simplest timestepper for differential equations. It defines the solution at the next time step, , in terms of the same at the previous time step, , using one-step forward-time finite different approximation of the time derivative, which for the system (9) gives
The time discretization step is presumed here the same for all steps of a simulation.
Ii-B2 Matrix Rush-Larsen
The proposed Matrix Rush-Larsen method (MRL) assumes that the matrix changes only slightly during one time step and therefore can be approximated by a constant. The solution of (9) can then be written in terms of the matrix exponential,
We assume that the matrix is diagonalizable, i.e. can be represented in the form , where matrix is composed of the eigenvectors concatenated as column vectors, and matrix contains eigenvalues placed on the corresponding places on the diagonal. A sufficient condition of diagonalizability of a matrix is that all its eigenvalues are distinct, and this is the generic situation; but we of course check that it actually takes place in every case. Then the matrix exponential is calculated as
where the exponential of the diagonal matrix is obtained by exponentiation of its diagonal elements.
As the numerical solution of the eigenvalue problem is computationally expensive, we precompute the matrices , and for a fine grid of physiological potentials, , , , , , (all in mV) before compile time and save them in a file.
At start time, the eigenvalue and eigenvector matrices are loaded from the file and we precompute, for used in the particular simulation, the transition matrices
for all . At the run time, the solver simply refers to the tabulated transition matrix ,
where is the tabulated transmembrane potential that is the nearest to .
Along with the code, we provide precomputed files for a voltage step size of (size of 4.85 MB), that are sufficient to obtain accurate results. The tables with of 48.5 MB size, used for the simulations presented, are available from the authors upon request.
The method of tabulation (tab.) can be applied to all the presented numerical methods. However, its benefit is most essential in the MRL method, because matrix exponentiation is computationally expensive. The accuracy of the tabulation is dependent on the voltage step (here 0.01 mV) which is a matter of choice depending on memory availability and allowable pre-compile and start-time computation time.
Ii-B3 Hybrid Operator Splitting
The MRL achieves the purpose in principle but is relatively costly as multiplication by a dense matrix is required at each time step. On the other hand, it did not at all exploit the specific structure of the TR, illustrated by Fig. 3, that is, that the matrix is sparse and some TRs are much faster than others for some voltage ranges. Hence we propose a hybrid operator splitting method (HOS), which combines FE and MRL, and exploits the asymptotic structure of the TRs. In this method, we set
as described in Fig. 3: contains only TRs that are fast at high values of (, , , , and ); contains only TRs that are fast at low values of (, , , and ); and contains only uniformly slow TRs (, , , , , , , , , and ). Explicit expressions for are given in the Supplement.
Every timestep is then done in three substeps,
In our case the matrix exponentials in the two fast subsystems (17) and (18) are found analytically, through solving the corresponding ODE systems. This is possible because some of the equations corresponding to the matrices and are coupled in a specific manner and can be solved one by one where solution of one equation is substituted in the next etc. The full expressions and the method of derivation are given in the Supplement; here we present the solution for state in the equation (17) as an example:
The slow subsystem (19) uses FE, and since it contains only uniformly slow TRs, it can tolerate large time-steps, allowed by other components of the cell model, without loss of stability.
Ii-C A priori error estimates
Estimates by standard methods (e.g. [14, Chapter 5], see details in the Supplement), show that all three numerical schemes have local truncation errors of the second order, i.e. , although the coefficients vary: for FE we have , for MRL we have , and for HOS it is composed of contributions of the three substeps plus the error due to operator splitting, , where , and . So comparison of MRL and HOS with FE depends on the solution, but in any case accuracy of HOS it contingent on , and not being large at the same time, to ensure relative smallness of .
Most of the algorithms described here were implemented in C language in double precision floating point arithmetics and compiled using GNU Compiler Collection (GCC) (version 4.7.2). The exception is computation of eigenvalues and eigenvector tables, which was done using mathematical software Sage  (version 5.9). Simulation were performed on Intel Core i5-3470 CPU with the clock frequency 3.20GHz under GNU/Linux operating system (distribution Fedora 18).
Figure 4 shows the detail of the first millisecond of simulated cardiac excitation, the onset of an action potential (AP). The Markov chain model was solved using the three suggested integration methods: forward Euler (FE), matrix Rush-Larsen (MRL), and hybrid operator splitting (HOS), as described in the methods section. The model was solved with time step , , and , except for FE, which was also solved for , to be used as a reference, but not for , due to instability.
The model excitable cell was initially at the resting state, and at the time , an AP was initialized by instantaneous injection of potassium ions, raising the membrane potential to . The initial conditions of the states of the MC model are specified in Tab. I and the initial states of the remaining variables of the model can be found in the supplementary material. Before the initiation, more than 90% of the channels reside in the states and , which require at least three transitions to get to the open state. After the initiation, the channels start to transit rapidly towards the open state and then to the state . Within about almost all channels reside in the state . Then, the channels slowly transit to the state , where they stay until the resting potential is recovered. The states and are similar to the situation when the inactivation gate is closed in the gate model. The states and have less than 10% occupancy during all the stages of the action potential.The plot of is omitted, as this variable changes very little during the time interval shown.
The results for the time step are consistent in all panels. The FE is still stable at time step , however, compared to MRL and HOS, the FE solution is less accurate, resulting in a higher peak and faster decay of both the open state occupancy and the resulting current.
Comparison of the solutions for with the reference , obtained by FE with , is shown on Fig. 4 (first column, fourth row). We see that MRL and HOS approximate better than FE at the same time steps. This is consistent with results of evaluation of the error estimates over the AP solution: we have , , , with , all in , and , (see the Supplement). This suggests that exponential integrators can be useful, for their accuracy, even when instability is not a concern, say in systems with slower dynamics, such as .
At longer time steps, FE is unstable (Fig. 1 illustrates a mild case of the instability), while MRL and HOS continue to provide stable solutions. At , the peak of the most important component of , the occupancy of the open state , is slightly lower than at shorter timesteps. On the other hand, the decrease of the peak of the total current in these two methods (HOS, MRL at ) is relatively small compared to the decrease of the open state occupancy. Also, the decay of the current in the MRL is slower than in the other cases. Note that the lead of the APs onsets at against smaller time steps is comparable to the value of .
Approximation of the whole APs rather than just their onsets is illustrated in Fig. 5.
Further increase of the time steps (not shown) in MRL and HOS gives significant errors in the AP, e.g. at there is a overshoot. Stability persists for much longer: for HOS the solution becomes unphysical (a negative concentration) at about without loss of stability, and for MRL an instability occurs at about , although the solution is then also very different from the true AP.
Table III illustrates the efficiency of the three methods at three different time steps . This was done by measuring time taken by simulations consisting of 100 pulses with a cycle length (CL) of without any output. The pulses were initialized by an instantaneous injection of potassium ions of a sufficient amount to set the membrane potential to . The table shows times taken by the whole cell model (“Total”) and by the Markov Chain model computations (“”). The times shown are median values from six separate simulations in each case to minimize the effect of other processes running on the computer.
At , FE is the most efficient method. Computation of accounts for 21.8% and 31.5% of the overall computation cost in FE and HOS respectively. Tabulation allows reduction of the computation cost of by 49.1% and 65.3% in FE and HOS. MRL was used only with tabulation using the precomputed eigenvalues and eigenvectors matrices and the computational cost at the is comparable with FE. These proportions are consistent with the results at time step and for MRL and HOS. So, at the same time step, the computational costs of the proposed methods are slightly higher, but the accuracies are somewhat better, compared to FE. The most important benefit of HOS and MRL is, however, the possibility of using larger time steps.
Both proposed methods maintain stability at larger time steps, and improve the accuracy of the solution at the same time step, compared to the explicit ODE solver (FE). When tabulated, those methods are comparable to FE in computational cost. As expected, using larger time steps results in reduction of computational cost.
MRL method extends the popular RL method, developed for gate models, to Markov chain models. MRL is more universal than HOS, and may be made “automatic”. The only restriction of our implementation is the assumption of diagonalizability of matrix for all voltages. If in another model this happens not to be the case, then some more sophisticated approach would be needed. If non-diagonalizability is a regular feature, say due to identical definitions of some of the TRs, then a Jordan form can be used instead; if it only happens at selected voltages, then interpolation of matrices may be sufficient.
HOS method depends on the possibility to split the transition rates to multiple (three in our case) sub-systems according to their speeds, and solve each of the subsystems on its own. Our solution benefits from the possibility of solving the fast subsystems analytically. Implementing the analytical solution results in even better speed-up as the resulting timestepping matrices are sparse. However, the possibility of a suitable analytical solution is not guaranteed for a general MC model. In this case, the fast time subsystems can be solved using diagonalization like in the MRL method, which might require additional computational time.
Finally we comment on the order of approximation. In this paper we considered first order schemes, and they are most popular in practice. However, the approximation order can be improved by using more sophisticated methods, both for the whole cell model (say using Runge-Kutta approach) and for the exponential solvers. For the original Rush-Larsen scheme, higher-order variants have been proposed and tested [16, 17], and the same ideas can be extended to the matrix case as well. Naturally, HOS method may then need to involve a more sophisticated operator splitting method to correspond.
Another appealing direction for further research is application of the proposed methods to other important MC models. MRL is straightforward for any MC where TRs depend only on one variable, otherwise tabulation will be a bit more problematic. The success of the HOS approach will depend on the asymptotic properties of the TRs.
During this study, TS was funded by University of Exeter PhD Studentship and VNB was partly supported by EPSRC grant EP/I029664/1. We are grateful to C.E. Clancy for the permission to use the original authors code for this study and S. Sherwin for encouraging discussions.
-  A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve,” J. Physiol, vol. 117, no. 4, pp. 500–544, 1952.
-  R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed. Boston: Brooks/Cole, 2011.
-  S. Rush and H. Larsen, “A practical algorithm for solving dynamic membrane equations,” IEEE Trans. Biomed. Eng., vol. 25, no. 4, pp. 389–392, 1978.
-  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, no. 10, pp. 2546–2548, 2009.
-  M. Perego and A. Veneziani, “An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations,” Electronic Transactions on Numerical Analysis, vol. 35, pp. 234–256, 2009.
-  M. E. Marsh, S. T. Ziaratgahi, and R. J. Spiteri, “The secrets to the success of the Rush-Larsen method and its generalizations,” IEEE Trans. Biomed. Eng., vol. 59, no. 9, pp. 2506–2515, 2012.
-  A. Reidman and K. Trivedi, “Numerical transient analysis of Markov models,” Comput. Opns Res., vol. 15, no. 1, pp. 19–36, 1988.
-  C. E. Clancy and Y. Rudy, “Na channel mutation that causes both Brugada and long-QT syndrome phenotypes: a simulation study of mechanism,” Circulation, vol. 105, no. 10, pp. 1208–1213, 2002.
-  C. H. Luo and Y. Rudy, “A dynamic-model of the cardiac ventricular action-potential. 1. Simulations of ionic currents and concentration changes,” Circulation Research, vol. 74, no. 6, pp. 1071–1096, 1994.
-  J. L. Zeng, K. R. Laurita, D. S. Rosenbaum, and Y. Rudy, “Two components of the delayed rectifier K current in ventricular myocytes of the guinea-pig type. theoretical formulation and their role in repolarization,” Circulation Research, vol. 77, no. 1, pp. 140–152, 1995.
-  P. C. Viswanathan, R. M. Shaw, and Y. Rudy, “Effects of and heterogeneity on action potential duration and its rate dependence. a simulation study,” Circulation, vol. 99, no. 18, May 1999.
-  W. Stein et al., Sage Mathematics Software (Version 5.9), The Sage Development Team, 2013, http://www.sagemath.org.
Exponential integrators for a Markov chain model
of the fast sodium channel of cardiocytes
Tomáš Starý, Vadim N. Biktashev
The numeration of equations in this document continues from the main text, and the literature references are to the literature list in the main text, repeated in the end of this document for the reader’s convenience.
Appendix A Cell Model Definition
This section contains the definition of the model according to the authors’ code . The format of equations and subsections aims to correspond to the papers where those equations were published to facilitate a straightforward comparison. The known differences with the papers are marked by the sign:. Voltages are measured in mV, time in ms and concentrations in mmol/L.The membrane currents are adjusted to the specific membrane capacitance  and are measured in .
Standard ionic concentrations
which differs from  where ; .
Initial Values of Variables
which differs from  where ; , and no initial values were given for and .
- pump :
which is identical to 
, the Slow Component of the Delayed Rectifier Current
The definition of in equation (53) is multiplied by 0.615 “to simulate the intramural heterogeneity”, which is slightly different from the factor used in Viswanathan et al. (1999)  to simulate epicardial cell.
Otherwise the definition is identical to .
, the Fast Component of the Delayed Rectifier Current
which is identical to . The original notation for was ; we use for the gas constant in this section, and for one of the state occupancies of the Markov Chain model elsewhere in the rest of the paper.
which is identical to .