On the Computational Complexity of the Bidomain and the Monodomain Models of Electrophysiology

Abstract

The bidomain model, coupled with accurate models of cell membrane kinetics, is generally believed to provide a reasonable basis for numerical simulations of cardiac electrophysiology. Because of changes occurring in very short time intervals and over small spatial domains, discretized versions of these models must be solved on fine computational grids, and small time-steps must be applied. This leads to huge computational challenges that have been addressed by several authors. One popular way of reducing the CPU demands is to approximate the bidomain model by the monodomain model, and thus reducing a two by two set of partial differential equations to one scalar partial differential equation; both of which are coupled to a set of ordinary differential equations modeling the cell membrane kinetics. A reduction in CPU time of two orders of magnitude has been reported. It is the purpose of the present paper to provide arguments that such a reduction is not present when order-optimal numerical methods are applied. Theoretical considerations and numerical experiments indicate that the reduction factor of the CPU requirements from bidomain to monodomain computations, using order-optimal methods, typically is about 10 for simple cell models and less than two for more complex cell models.

Appendices

APPENDIX A: THE FITZHUGH–NAGUMO MODEL

We have used a modified FitzHugh–Nagumo model, given by

$$\begin{array}{l} \displaystyle I_{{\rm ion}} = k(v - v_{{\rm rest}} )(v - v_{{\rm th}} )(v - v_{{\rm peak}} ) - k(v - v_{{\rm rest}} )w, \\ \displaystyle \frac{{dw}}{{dt}} = l(v - v_{{\rm rest}} - bw). \\\end{array}$$

The values of the parameters are given in Table 6.

TABLE 6. Values of the parameters in the modified FitzHugh–Nagumo model.

APPENDIX B: THE BEELER–REUTER MODEL

The Beeler–Reuter model describes four different ionic currents, which are controlled by six gating variables. The model also describes the intracellular calcium concentration, resulting in a total of eight ODEs:

$$I_{{\rm ion}} = I_{{\rm Na}} + I_{\rm K} + I_x + I_s,$$

(27)

$$\frac{{dc}}{{dt}} = 0.07(1 - c) - I_s,\vspace*{-12pt}$$

(28)

$$\frac{{dg}}{{dt}} = \alpha _g (1 - g) - \beta _g g,$$

(29)

with g=m, h, j, d, f, x.

For notational convenience the intracellular calcium concentration has been scaled, c=10⁷[Ca] _i . The ionic currents are given by

$$\begin{array}{l} \displaystyle I_{{\rm Na}} = (4m^3 hj + 0.033)(v - 50), \\ \displaystyle I_{\rm K} = 1.4\frac{{\exp (0.04(v + 85)) - 1}}{{\exp (0.08(v + 53)) + \exp (0.04(v + 53))}}, \\ \displaystyle I_x = 0.8x\frac{{\exp (0.04(v + 77)) - 1}}{{\exp (0.04(v + 35))}}, \\ \displaystyle I_s = 0.09fd(v + 66.18 + 13.0287\ln [{\rm Ca}]_i ).\end{array}$$

The rate functions α _g , β _g in (29) are all of the form

$$\alpha _g = \frac{{C_1 \exp (C_2 (v - v_0 )) + C_3 (v - v_0 )}}{{1 + C_4 \exp (C_5 (v - v_0 ))}},$$

(30)

with the constants specified in Table 7.

TABLE 7. The constants in the rate functions of the Beeler–Reuter model.

APPENDIX C: THE MODEL BY WINSLOW et al.

The model by Winslow et al. is substantially more complex than the Beeler–Reuter model. The model consists of 33 ODEs, and the total ionic current is defined as a sum of 13 different ionic currents. The specification of the complex ionic current and rate functions goes over several pages, and is not included here. A complete specification of the model can be found in the original publication,23 and a brief summary is also given in Sundnes et al. 17