Incorporating inductances in tissue-scale models of cardiac electrophysiology

Rossi, Simone; Griffith, Boyce E.

doi:10.1063/1.5000706

In standard models of cardiac electrophysiology, including the bidomain and monodomain models, local perturbations can propagate at infinite speed. We address this unrealistic property by developing a hyperbolic bidomain model that is based on a generalization of Ohm’s law with a Cattaneo-type model for the fluxes. Further, we obtain a hyperbolic monodomain model in the case that the intracellular and extracellular conductivity tensors have the same anisotropy ratio. In one spatial dimension, the hyperbolic monodomain model is equivalent to a cable model that includes axial inductances, and the relaxation times of the Cattaneo fluxes are strictly related to these inductances. A purely linear analysis shows that the inductances are negligible, but models of cardiac electrophysiology are highly nonlinear, and linear predictions may not capture the fully nonlinear dynamics. In fact, contrary to the linear analysis, we show that for simple nonlinear ionic models, an increase in conduction velocity is obtained for small and moderate values of the relaxation time. A similar behavior is also demonstrated with biophysically detailed ionic models. Using the Fenton–Karma model along with a low-order finite element spatial discretization, we numerically analyze differences between the standard monodomain model and the hyperbolic monodomain model. In a simple benchmark test, we show that the propagation of the action potential is strongly influenced by the alignment of the fibers with respect to the mesh in both the parabolic and hyperbolic models when using relatively coarse spatial discretizations. Accurate predictions of the conduction velocity require computational mesh spacings on the order of a single cardiac cell. We also compare the two formulations in the case of spiral break up and atrial fibrillation in an anatomically detailed model of the left atrium, and we examine the effect of intracellular and extracellular inductances on the virtual electrode phenomenon.

Since its introduction by Hodgkin and Huxley,³⁷ the cable equation and its higher-dimensional generalizations⁸³ have been common models of electrical impulse propagation in excitable media, including neurons and muscle. The effects of inductances in these systems are considered to be relatively small, and so they are neglected in classical versions of these models. By omitting such terms, the standard equations of cardiac electrophysiology become parabolic, but, as in all parabolic equations, local perturbations can propagate at infinite speed. This unrealistic property has been addressed in models of neurons by Lieberstein and Mahrous,⁵³ and hyperbolic models including inductances have been proposed by Engelbrecht,²¹ Lieberstein,⁵² and Engelbrecht et al.²² These models are also supported by the fact that neurons, skeletal muscle cells, and cardiomyocytes show typical resonance effects due to inductances, as demonstrated by the studies of Clapham and Defelice¹⁵ and Koch.⁴⁶ The common conclusion that inductances are negligible, which is based on the linear analyses for neurons by Scott^74,75 and the numerical results of Kaplan and Trujillo,⁴¹ may not be valid in the complex arrangement of cardiac tissue, where the inhomogeneities together with highly nonlinear reactions can lead to reentrant waves and chaotic behavior. We thus derive a hyperbolic model for cardiac electrophysiology, and we compare the solutions with parabolic models in several cases, including simple excitation patterns as well as for spiral waves, atrial fibrillation, and the virtual electrode phenomenon.

The heart is a complex organ with a highly heterogeneous structure. Muscle cells are embedded in an extracellular compartment that includes many components, including capillaries, connective tissue, and collagen. The structural arrangement of the tissue is known to influence electrical impulse propagation.^12,78,84 The standard models of electrophysiology neglect all these complexities, and the resulting equations have the unrealistic feature that local perturbations can propagate infinitely fast. These models are fundamentally based on the cable equation. The studies of Lieberstein⁵² and Lieberstein and Mahrous⁵³ were the first to suggest that the cable model for neurons should contain inductances because of the three-dimensional nature of the axon. Several authors, including Kaplan and Trujillo,⁴¹ Scott,⁷⁴ and Engelbrecht,²¹ have investigated this hypothesis for nerves. Their conclusions were that inductances in neurons are of the order of few μH and therefore are negligible, following a linear analysis of a one-dimensional nerve discussed by Scott.⁷⁵ Several years later, further studies by DeHaan and DeFelice,²⁰ Clapham and Defelice,¹⁵ and Koch⁴⁶ showed that the cell membranes of excitable tissue exhibit self-inductance. In particular, Clapham and Defelice¹⁵ showed that the impedance of the embryonic heart cell membrane resonates at a frequency around 1 Hz, thereby enhancing homogeneity of the voltage. To our knowledge, there has been no experimental study that addresses the role of inductances for reentrant waves. In this chaotic scenario, inductances may have an important role in the system dynamics.

Excitable tissues have a characteristic speed of transmission that limits the velocity at which signals can propagate. Any signal, including action potential wavefronts, cannot propagate faster than this characteristic speed. On the other hand, it is important to notice that propagation of the action potential is related to the nonlinear dynamics of the system, and not to the wave-like behavior that may be induced by inductances. One way to transform the cable equation into a hyperbolic equation is to simply add an “inertial” term proportional to the persistence time of the diffusive process, as done by Zemskov et al.⁸⁷ Unfortunately, this type of model cannot be derived from reasonable physical assumptions. However, as we show here, it is possible to perform a phenomenological derivation of a hyperbolic reaction-diffusion model by using the Cattaneo model for the fluxes. This formulation was originally introduced by Cattaneo¹⁴ to eliminate the anomalies found using Fourier’s law in the heat equation, and this model has subsequently been used in a wide range of applications, including forest fire models,⁵⁷ chemical systems,^29,51 thermal combustion,²⁶ and the spread of viral infections.²⁷ Cattaneo-type models for the fluxes have been derived in several ways, ranging from phenomenological and thermodynamical derivations to isotropic and anisotropic random walks with reactions.^40,56 The use of Cattaneo-type fluxes in a monodomain model of cardiac electrophysiology leads to two additional terms in the equations proportional to the characteristic relaxation time of the medium: the second derivative in time of the potential, which is associated with “inertia,” and the time derivative of the ionic currents. Even if the relaxation time is small, the rapid variation of the ionic currents can impart a contribution that may not be negligible. This is particularly relevant near the front of the wave, where fast currents give rise to the upstroke of the action potential.

Verification and validation remain major challenges in computational electrophysiology. Krishnamoorthi et al.⁴⁸ propose that the “wave speed should not be sensitive to choices of the numerical solution protocol, such as mesh density, numerical integration scheme, etc.” Although this is a fundamentally desirable criterion, it is also nearly impossible to achieve in practice. What does seem reasonable is to require that that the error in the wave speed be “small,” but how small the error must be taken clearly depends on the case under consideration. For example, for a single heart beat, an error of 5% might represent a reasonable approximation. On the other hand, models of cardiac fibrillation may require a much smaller error. In fact, in this case, a 5% error in the conduction velocity can determine whether reentrant waves form, or when and how they break up. Although this problem has been discussed in detail in prior work,⁶⁶ the common belief within the field seems to be that spatial discretizations on the order of 200 μm are sufficient to capture the conduction velocities of the propagating fronts. This estimate seems to hold for isotropic propagation at normal coupling strengths, but in many important cases, the most relevant propagation of the electrical signal is transverse to the alignment of the cardiac cells, where the conductivities are typically 8 times smaller than in the longitudinal direction.³³ As shown by Quarteroni et al.,⁶⁸ the mesh sizes needed to resolve transverse propagation are actually closer to 25 μm. Such a small mesh spacing requires the use of large-scale simulations and highly efficient codes. Further, the need to use such high spatial resolutions challenges the fundamental idea of using a continuum model for the description of cardiac electrophysiology, and multiscale models have been proposed.³² This paper shows that the transverse conduction velocities are sensitive to the grid size and to the mesh orientation for both regular and hyperbolic versions of the model, and that the hyperbolic model has similar mesh size requirements as the standard model.

I. PHENOMENOLOGICAL DERIVATION OF THE HYPERBOLIC BIDOMAIN EQUATIONS

In the 1970s, Tung⁸³ formulated a bidomain model of the propagation of the action potential in cardiac muscle. This tissue-scale model considers the myocardium to be composed of continuous intracellular and extracellular compartments, coupled via a continuous cellular membrane. Although the bidomain equations can be derived from a model that accounts for the tissue microarchitecture,^33,42,60 the bidomain model is a fundamentally homogenized description of excitation propagation that neglects the details of this microarchitecture. Instead, the bidomain equations describe the dynamics of a local average of the voltages in the intracellular and extracellular compartments over a control volume. One of the assumptions required by the homogenization procedure is that the control volume is large compared to the scale of the cellular microarchitecture but small compared to any other important spatial scale of the system, such as the width of the action potential wavefront. Although the validity of this model has been questioned, for example by Bueno-Orovio et al.,¹³ this approach has been extremely successful, and at present, most simulations of cardiac electrophysiology use such models. For a detailed review of the bidomain model and other models of electrophysiology, we refer to Griffith and Peskin.³¹ Here, we assume that the homogenization assumptions hold, and we derive the hyperbolic bidomain model phenomenologically, starting from charge conservation in quasistatic conditions.

In the hyperbolic bidomain model, as in the standard bidomain model, the extracellular and intracellular compartments are characterized by the anisotropic conductivity tensors, respectively, $σ_{e}$ and $σ_{i}$ ⁠. Introducing a local orthonormal basis, ${f_{0}, s_{0}, n_{0}}$ ⁠, we assume the conductivity tensors can be written as $σ_{j} = σ_{j}^{f} f_{0} \otimes f_{0} + σ_{j}^{s} s_{0} \otimes s_{0} + σ_{j}^{n} n_{0} \otimes n_{0},$ for $j = e, i$ ⁠. Typical values of the conductivity coefficients range from $10^{- 2}$ mS/cm in the cross-fiber direction to 1 mS/cm in the fiber direction (for example, Colli Franzone et al.¹⁷ use $σ^{f} = 1.2$ mS/cm, $σ^{s} = 0.346$ mS/cm, and $σ^{n} = 0.0435$ mS/cm). We model the current density fluxes using a Cattaneo-type equation, so that

τ_{e} \frac{\partial J_{e}}{\partial t} + J_{e} = - σ_{e} \nabla V_{e},

(1)

τ_{i} \frac{\partial J_{i}}{\partial t} + J_{i} = - σ_{i} \nabla V_{i},

(2)

where $V_{i}$ (⁠ $V_{e}$ ⁠), $J_{i}$ (⁠ $J_{e}$ ⁠), and $τ_{i}$ (⁠ $τ_{e}$ ⁠) are the intracellular (extracellular) potential, the intracellular (extracellular) flux, and the intracellular (extracellular) relaxation time, respectively. As shown by Jou et al.,⁴⁰ relations (1) and (2) can be also derived from the simplest model of ionic conduction in a dilute system. The derivation of these evolution laws for the fluxes from the generalized Gibbs equation^28,40 shows that they are consistent with the theory of extended irreversible thermodynamics. Alternatively, Eqs. (1) and (2) can be interpreted as arising from a circuit model that includes inductances, such as the one depicted in Fig. 1. The derivation of the cable equation for the circuit model shown in Fig. 1 is performed in Appendix A. In particular, we show in Eq. (A15) the relationship between inductance and the relaxation time.

FIG. 1.

View large Download slide

Schematic diagram of a discretized cable including inductance effects, with isopotential circuit elements of length $Δ x$ ⁠. The inductances in the cable give a Cattaneo-type model of the fluxes of the form (1) and (2).

To derive the higher-dimensional model equations, we begin by taking the divergence of the fluxes (1) and (2), so that

τ_{e} \frac{\partial}{\partial t} \nabla \cdot J_{e} + \nabla \cdot J_{e} = - \nabla \cdot σ_{e} \nabla V_{e},

(3)

τ_{i} \frac{\partial}{\partial t} \nabla \cdot J_{i} + \nabla \cdot J_{i} = - \nabla \cdot σ_{i} \nabla V_{i} .

(4)

As in similar derivations of the bidomain model, we impose a quasistatic form of charge conservation,⁷² yielding

\nabla \cdot (J_{i} + J_{e}) = 0.

(5)

Additionally, the current leaving each compartment needs to enter the other, so that

- \nabla \cdot J_{i} = I_{t} = \nabla \cdot J_{e},

(6)

where $I_{t} = χ (C_{m} \frac{\partial V}{\partial t} + I_{ion})$ is the usual transmembrane current density, with χ being the membrane area per unit volume of tissue, $C_{m}$ the membrane capacitance, and $I_{ion}$ the transmembrane ionic current. Using Eq. (6) in (3) and (4), we obtain the system of equations

τ_{e} \frac{\partial I_{t}}{\partial t} + I_{t} = - \nabla \cdot σ_{e} \nabla V_{e},

(7)

- τ_{i} \frac{\partial I_{t}}{\partial t} - I_{t} = - \nabla \cdot σ_{i} \nabla V_{i} .

(8)

Defining the transmembrane potential $V = V_{i} - V_{e}$ to eliminate $V_{i}$ from the equations yields

τ_{e} \frac{\partial I_{t}}{\partial t} + I_{t} = - \nabla \cdot σ_{e} \nabla V_{e},

(9)

τ_{i} \frac{\partial I_{t}}{\partial t} + I_{t} = \nabla \cdot σ_{i} \nabla V + \nabla \cdot σ_{i} \nabla V_{e} .

(10)

Expanding the transmembrane currents, we finally obtain the hyperbolic bidomain model

\begin{matrix} τ_{e} C_{m} \frac{\partial^{2} V}{\partial t^{2}} + C_{m} \frac{\partial V}{\partial t} + \nabla \cdot D_{e} \nabla V_{e} = - I_{ion} - τ_{e} \frac{\partial I_{ion}}{\partial t}, \end{matrix}

(11)

\begin{matrix} τ_{i} C_{m} \frac{\partial^{2} V}{\partial t^{2}} + C_{m} \frac{\partial V}{\partial t} - \nabla \cdot D_{i} \nabla V - \nabla \cdot D_{i} \nabla V_{e} \\ = - I_{ion} - τ_{i} \frac{\partial I_{ion}}{\partial t}, \end{matrix}

(12)

where we set $D_{i} = σ_{i} / χ$ and $D_{e} = σ_{e} / χ$ ⁠. Alternatively, the model equations can be written as

\begin{matrix} τ_{i} C_{m} \frac{\partial^{2} V}{\partial t^{2}} + C_{m} \frac{\partial V}{\partial t} - \nabla \cdot D_{i} \nabla V - \nabla \cdot D_{i} \nabla V_{e} \\ = - I_{ion} - τ_{i} \frac{\partial I_{ion}}{\partial t}, \end{matrix}

(13)

\begin{matrix} (τ_{e} - τ_{i}) C_{m} \frac{\partial^{2} V}{\partial t^{2}} + \nabla \cdot D_{i} \nabla V + \nabla \cdot (D_{e} + D_{i}) \nabla V_{e} \\ = (τ_{i} - τ_{e}) \frac{\partial I_{ion}}{\partial t} . \end{matrix}

(14)

Notice that if $τ_{i} = τ_{e} = τ$ ⁠, then these equations reduce to a hyperbolic-elliptic system that is similar to the parabolic-elliptic form of the standard bidomain model

\begin{matrix} τ C_{m} \frac{\partial^{2} V}{\partial t^{2}} + C_{m} \frac{\partial V}{\partial t} - \nabla \cdot D_{i} \nabla V - \nabla \cdot D_{i} \nabla V_{e} = - I_{ion} - τ \frac{\partial I_{ion}}{\partial t}, \end{matrix}

(15)

\nabla \cdot D_{i} \nabla V + \nabla \cdot (D_{e} + D_{i}) \nabla V_{e} = 0.

(16)

If we further take τ = 0, we retrieve the usual bidomain model in its parabolic-elliptic form

C_{m} \frac{\partial V}{\partial t} - \nabla \cdot D_{i} \nabla V - \nabla \cdot D_{i} \nabla V_{e} = - I_{ion},

(17)

\nabla \cdot D_{i} \nabla V + \nabla \cdot (D_{e} + D_{i}) \nabla V_{e} = 0.

(18)

II. REDUCTION TO THE HYPERBOLIC MONODOMAIN

The hyperbolic bidomain model can be simplified by assuming that the extracellular and intracellular compartments have the same anisotropy ratios, so that $D = D_{e} = λ D_{i}$ ⁠. If we make this assumption, we obtain from (14)

\begin{matrix} - \nabla \cdot D \nabla V_{e} = \frac{1}{λ + 1} \nabla \cdot D \nabla V + \frac{λ}{λ + 1} (τ_{e} - τ_{i}) C_{m} \frac{\partial^{2} V}{\partial t^{2}} - \frac{λ}{λ + 1} (τ_{i} - τ_{e}) \frac{\partial I_{ion}}{\partial t} . \end{matrix}

(19)

Substituting in (13), we find

\begin{matrix} [τ_{i} + \frac{λ}{λ + 1} (τ_{e} - τ_{i})] C_{m} \frac{\partial^{2} V}{\partial t^{2}} + C_{m} \frac{\partial V}{\partial t} - \frac{λ}{λ + 1} \nabla \cdot D \nabla V \\ = - I_{ion} - [τ_{i} + \frac{λ}{λ + 1} (τ_{e} - τ_{i})] \frac{\partial I_{ion}}{\partial t} . \end{matrix}

(20)

Defining $τ = τ_{i} + λ (τ_{e} - τ_{i}) / (λ + 1)$ ⁠, we obtain the hyperbolic monodomain model

τ C_{m} \frac{\partial^{2} V}{\partial t^{2}} + C_{m} \frac{\partial V}{\partial t} - \nabla \cdot D \nabla V = - I_{ion} - τ \frac{\partial I_{ion}}{\partial t},

(21)

where here we have absorbed the term $λ / (λ + 1)$ into D. The relaxation time τ of the monodomain model is always positive, and it is zero only if both $τ_{i}$ and $τ_{e}$ are zero. In fact, if $τ_{e} = 0$ ⁠, then $τ = τ_{i} / (λ + 1)$ ⁠, and if $τ_{i} = 0$ ⁠, then $τ = λ τ_{e} / (λ + 1)$ ⁠. However, if $τ_{i} = τ_{e}$ ⁠, then $τ = τ_{i} = τ_{e}$ ⁠.

Introducing a new variable $Q = \frac{\partial V}{\partial t}$ ⁠, we transform the hyperbolic monodomain equations into the first-order system

\frac{\partial V}{\partial t} = Q,

(22)

τ C_{m} \frac{\partial Q}{\partial t} + C_{m} Q - \nabla \cdot D \nabla V = - I_{ion} - τ \frac{\partial I_{ion}}{\partial t} .

(23)

Equations (22) and (23) are usually supplemented with insulation boundary conditions, such that

\nabla V \cdot N = 0

(24)

on the boundary of the domain, where the vector N is the normal to the boundary.

Following Stan et al.,⁸⁰ we say, that a solution of the hyperbolic monodomain equations has a finite propagation speed if, given compactly supported initial conditions for V at time t = 0, $V (\cdot, t)$ is also compactly supported for any t > 0. The compact support is taken with respect to the resting potential V₀. By contrast, a solution has infinite propagation speed if the initial data are compactly supported, but for any t > 0 and any R > 0, the set $M_{R, t} = {x : | | x | |_{2} \geq R, V (x, t) > V_{0}}$ has a positive measure. For any relaxation time $τ \neq 0$ ⁠, (21) is hyperbolic, and it thereby has the property of finite propagation speed. Setting τ = 0, the equations become parabolic and have solutions with infinite propagation speeds. Specifically, in the standard parabolic model, local perturbations in V, even those that do not generate a propagating front, will travel at infinite speed.

III. IONIC MODELS

Transmembrane ionic fluxes through ion channels, pumps, and exchangers are responsible for the cardiac action potential. The action potential is initiated by a fast inward sodium current that depolarizes the cellular membrane. After depolarization phase, slow inward currents (primarily calcium currents) and slow outward currents (primarily potassium currents) approximately balance each other, prolonging the action potential and creating a plateau phase. Ultimately, the slow outward currents bring the transmembrane potential difference back to its resting value of approximately $- 80 mV$ ⁠. The bidomain and monodomain models must be completed by specifying the form of $I_{ion}$ ⁠, which accounts for these transmembrane currents. The states of the transmembrane ion channels are described by a collection of variables, w, associated with the ionic model, so that $I_{ion} = I_{ion} (V, w)$ ⁠. Typically, the state variable dynamics are determined by a spatially decoupled system of nonlinear ordinary differential equations

\frac{\partial w}{\partial t} = g (V, w),

(25)

where the form of g depends on the details of the particular ionic model. We consider the simplified piecewise-linear model of McKean,⁵⁵ the two-variable model of Aliev and Panfilov (AP),³ the three-variable model of Fenton and Karma,²⁵ and the biophysically detailed models of Tusscher and Panfilov⁸¹ for the ventricles (20 variables) and of Grandi et al.³⁰ for the atria (57 variables).

In the two-variable model of Aliev and Panfilov,³ $w = {r}$ ⁠, with

\frac{\partial r}{\partial t} = (ε + \frac{μ_{1} r}{μ_{2} + V}) (- r - k V (V - b - 1)),

(26)

and the ionic current is $I_{ion} = k V (V - α) (V - 1) + r V$ ⁠. The parameters used in this case are shown in Table I. In the three-variable model of Fenton and Karma,²⁵ $w = {v, w}$ ⁠, with

\frac{\partial v}{\partial t} = \frac{1}{τ_{v}^{-} (V)} (1 - p) (1 - v) - \frac{1}{τ_{v}^{+}} p v,

(27)

\frac{\partial w}{\partial t} = \frac{1}{τ_{w}^{-}} (1 - p) (1 - w) - \frac{1}{τ_{w}^{+}} p w,

(28)

where $τ_{v}^{-} (V) = (1 - q) τ_{v 1}^{-} + q τ_{v 2}^{-}, p = H (V - V_{c}), q = H (V - V_{v})$ ⁠, and $H (\cdot)$ is the Heaviside function. The total ionic current is the sum of three currents, $I_{ion} = - I_{fi} (V, v) - I_{so} (V) - I_{si} (V, w)$ ⁠, with

I_{fi} (V, v) = - \frac{1}{τ_{d}} v p (V - V_{c}) (1 - V),

(29)

I_{so} (V) = \frac{1}{τ_{0}} V (1 - p) + \frac{1}{τ_{r}} p,

(30)

I_{si} (V, w) = - \frac{1}{2 τ_{si}} w (1 + \tanh (k (V - V_{c}^{si}))) .

(31)

The sets of parameters used in this paper for this model are shown in Table II. We refer to the original papers^30,81 for the statements of the equations and parameters of the biophysically detailed ionic models.

TABLE I.

Parameters for the Aliev–Panfilov ionic model (26).

$C_{m}$	$σ^{f}$	$σ^{s} = σ^{n}$	χ	k	b	μ₁	μ₂	ε
1	1	0.125	1	8	0.1	0.12	0.3	0.01

$C_{m}$	$σ^{f}$	$σ^{s} = σ^{n}$	χ	k	b	μ₁	μ₂	ε
1	1	0.125	1	8	0.1	0.12	0.3	0.01

TABLE II.

Sets of parameters for the Fenton–Karma ionic model (27)–(31) as stated by Fenton et al.²⁴

	$C_{m}$	$σ^{f}$	$σ^{s} = σ^{n}$	χ	$τ_{v}^{+}$	$τ_{v 1}^{-}$	$τ_{v 2}^{-}$	$τ_{w}^{+}$	$τ_{w}^{-}$	$τ_{d}$	τ₀	$τ_{r}$	$τ_{si}$	k	$V_{c}^{si}$	$V_{c}$	V_v
Parameter set 3	1	0.1	0.0125	1	3.33	19.6	1250	870	41	0.25	12.5	33.33	29.0	10	0.85	0.13	0.04
Parameter set 4	1	0.1	0.0125	1	3.33	15.6	5	350	80	0.407	9.0	34.0	26.5	15	0.45	0.15	0.04
Parameter set 5	1	0.1	0.0125	1	3.33	12.0	2	1000	100	0.362	5.0	33.33	29.0	15	0.7	0.13	0.04
Parameter set 6	1	0.1	0.0125	1	3.33	9.0	8	250	60	0.395	9.0	33.33	29.0	15	0.5	0.13	0.04

	$C_{m}$	$σ^{f}$	$σ^{s} = σ^{n}$	χ	$τ_{v}^{+}$	$τ_{v 1}^{-}$	$τ_{v 2}^{-}$	$τ_{w}^{+}$	$τ_{w}^{-}$	$τ_{d}$	τ₀	$τ_{r}$	$τ_{si}$	k	$V_{c}^{si}$	$V_{c}$	V_v
Parameter set 3	1	0.1	0.0125	1	3.33	19.6	1250	870	41	0.25	12.5	33.33	29.0	10	0.85	0.13	0.04
Parameter set 4	1	0.1	0.0125	1	3.33	15.6	5	350	80	0.407	9.0	34.0	26.5	15	0.45	0.15	0.04
Parameter set 5	1	0.1	0.0125	1	3.33	12.0	2	1000	100	0.362	5.0	33.33	29.0	15	0.7	0.13	0.04
Parameter set 6	1	0.1	0.0125	1	3.33	9.0	8	250	60	0.395	9.0	33.33	29.0	15	0.5	0.13	0.04

IV. NUMERICAL RESULTS

The hyperbolic monodomain and bidomain models are discretized using a low-order finite element scheme described in Appendix C that is implemented in the open-source parallel C++ code BeatIt (available at http://github.com/rossisimone/beatit), which is based on the libMesh finite element library⁴⁴ and relies on linear solvers provided by PETSc.^6–8 All the codes used for the following tests are contained in the online repository, and all tests can be replicated directly from those codes. The only exception is the atrial fibrillation test, which uses a patient-specific mesh that is not contained in the repository. One- and two-dimensional simulations were run using a Linux workstation with two Intel Xeon E5–2650 v3 processors (up to 40 threads) and 32 GB of memory. Three-dimensional simulations were run on the KillDevil Linux cluster at the University of North Carolina at Chapel Hill. We used Matlab⁵⁴ to visualize the one-dimensional results and Paraview² for the two- and three-dimensional simulations.

A. Comparison with an exact solution

We start by considering a simple piecewise-linear bistable model for the ionic currents, with

\begin{matrix} I_{ion} (V) = k V - k [V_{2} H (V - V_{1}) + V_{0} (1 - H (V - V_{1}))], \end{matrix}

(32)

where V₀ is the resting potential, V₁ is the threshold potential, V₂ is the depolarization potential, and $H (\cdot)$ is the Heaviside function. This model (32) reduces to the piecewise-linear model of McKean⁵⁵ after a simple dimensional analysis, after which the nondimensional ionic currents take the form

{\hat{I}}_{ion} (\hat{V}) = \hat{V} - H (\hat{V} - α) .

(33)

The nondimensional potential $\hat{V}$ is related to the dimensional potential by $V = (V_{2} - V_{0}) \hat{V} + V_{0}$ ⁠. In this model, $α = (V_{1} - V_{0}) / (V_{2} - V_{0})$ describes the excitability of the tissue. As we show in Appendix B, using model (33), it is possible to find the analytic solution of a propagating front for the nondimensional hyperbolic monodomain problem. In particular, we find that the front propagation speed in an unbounded domain is

c = \frac{(1 - 2 α)}{\sqrt{μ + (α - α^{2}) {(μ - 1)}^{2}}} < \sqrt{\frac{1}{μ}} = c_{s} .

(34)

The speed $c_{s} = \sqrt{1 / μ} = \sqrt{C_{m} / τ k}$ represents the maximum speed at which a perturbation can travel in the system. When the relaxation time τ goes to zero, the local perturbations can travel at infinite speed. The parameter $μ = τ k / C_{m}$ is a nondimensional number representing the ratio between the relaxation time and the characteristic time of the reactions and characterizing the effects of the inductances in the system (see Appendix B). The corresponding dimensional speed is

v = \sqrt{\frac{σ k}{χ C_{m}^{2}}} \frac{(1 - 2 α)}{\sqrt{μ + (α - α^{2}) {(μ - 1)}^{2}}} < \sqrt{\frac{σ}{χ τ C_{m}}} .

(35)

As a verification test, we consider the spatial interval $Ω = [0, 50]$ ⁠. An initial stimulus of amplitude 1 is applied at $x \in [24.5, 25.5]$ for the interval $t \in [0.03, 1.03]$ ⁠. The system of equations is solved using a mesh size h = 0.03125 and a time step size $Δ t = 0.003653.$ We show the nondimensional solutions $\hat{V}$ and $\hat{Q}$ at t = 10 in Fig. 2 (center and right). We register the activation time at a particular spatial location whenever $\hat{V}$ crosses a threshold of 0.9. The conduction velocities are measured by picking the distance between two points x₁ and $x_{2} \in Ω$ and dividing it by the time interval between the activation times at these two locations. Because (34) is obtained by assuming that Ω is the entire real line, we need to measure the conduction velocities far from the boundaries. For this reason, we choose $x_{1} = 30$ and $x_{2} = 32$ ⁠. The comparison between the exact conduction velocities defined by Equation (34) and those obtained by the simulations is shown in Fig. 2 (left).

FIG. 2.

View large Download slide

Conduction velocities for the piecewise-linear reaction model (33). (Left) comparison between the numerical and the exact wave speed (34). The black curve represents the characteristic propagation speed of the system, i.e., the limit speed at which local perturbations can travel. If the relaxation time τ = 0, then local perturbations can propagate at infinite speed. (Center) propagating front for different relaxation times at t = 10. (Right) time derivative of the nondimensional potential for different relaxation times at t = 10. All simulations were run using a mesh size h = 0.03125 and a time step size $Δ t = 0.003653.$ ⁠.

We also examine the effect of the relaxation time on the conduction velocities for three values of the excitability parameter, $α \in {0.1, 0.2, 0.3}$ ⁠. With this simple ionic model, the larger the relaxation time, the slower the wave speed. The limit speed $c_{s}$ at which local perturbations can travel is also shown in Fig. 2 (left). Thus, in this piecewise-linear model, the effect of inductances is to slow down the propagation of the front. By contrast, the following tests will show that with fully nonlinear models, inductances can also enhance propagation.

Using the analytic solution derived in Appendix B, we also perform a space-time convergence study. We consider the spatial interval $Ω = [- 25, 25]$ and set the initial conditions according to (B9) and (B10), assuming the front at time t = 0 is at x = 0 and $α = 0.1$ ⁠. The system is discretized in time using implicit-explicit (IMEX) Runge–Kutta (RK) time integrators.^5,10 In particular, we use the forward/backward Euler schemes for the first-order time integrator and the Heun/Crank–Nicholson schemes for the second-order time integrators. For more details on the time integrator, see Appendix C. To ensure the validity of the solutions (B9) and (B10) in the considered bounded domain, we run the simulation from t = 0 to t = 1, so that the front is still far from the boundaries. The time step size $Δ t$ is taken to be linearly proportional to the mesh size h, with a value of $Δ t = 0.05$ in the coarsest cases. The time step size was chosen to guarantee a large enough number of time iterations while satisfying the Courant–Friedrichs–Lewy (CFL) condition (see Appendix C). Figure 3 shows the errors for the potential V and its time derivative Q using the first-order (Fig. 3, left) and second-order (Fig. 3, center) time stepping schemes. Because the ionic currents and their derivatives in this case are not smooth functions, we do not expect to obtain full second-order convergence. On the other hand, whereas the convergence rates for the parabolic monodomain model are always first order in both V and Q, the hyperbolic monodomain model with the second-order time stepping scheme converges quadratically in V and linearly only in Q. The simulations were run with and without regularization of the Heaviside and Dirac-δ functions. The errors reported in Fig. 3 correspond to the most accurate simulations. We show in Fig. 3 (right) that the slow convergence results from the non-smoothness of the ionic currents: we compare the convergence of the piecewise-linear reaction model with the monodomain model with the cubic reaction term

{\hat{I}}_{ion} (\hat{V}) = \hat{V} (\hat{V} - 1) (\hat{V} - α) .

The exact solution for the cubic reaction in the parabolic monodomain has been determined previously.^42,66 We are not aware of the existence of an analytical solution for the hyperbolic monodomain with a cubic reaction term. We show in Fig. 3 that the cubic model converges quadratically (using a second-order time stepping scheme) whereas the piecewise-linear model converges linearly. To conclude, at least in these tests, the discretization errors of the parabolic mondomain model are larger than the discretization errors in the hyperbolic monodomain model, indicating a better accuracy for the hyperbolic model.

FIG. 3.

View large Download slide

Convergence study for the piecewise-linear reaction model (33) for the variables V and Q. (Left) comparison of the rate of convergence for the hyperbolic and parabolic monodomain models using a first-order implicit-explicit (IMEX) Runge–Kutta (RK) time integrator. (Center) comparison of the rate of convergence for the hyperbolic and parabolic monodomain models using a second-order IMEX-RK time integrator. (Right) comparison of the rate of convergence between the piecewise-linear reaction model and a cubic reaction model for the parabolic monodomain model using a second-order IMEX-RK time integrator.

B. Effect of the relaxation time on the conduction velocity for different ionic models

We now analyze how the inductance terms influence the conduction velocities of cardiac action potentials. We consider four ionic models: the two-variable model of Aliev and Panfilov,³ the three-variable model of Fenton and Karma,²⁵ the ventricular M-cell model of Tusscher and Panfilov⁸¹ (20 variables), and the atrial model of Grandi et al.³⁰ (57 variables). For the two-variable Aliev–Panfilov model, we use the parameters of Nash and Panfilov,⁵⁹ as reported in Table I. We consider two values for the excitability parameter, $α = 0.1$ and $α = 0.2$ ⁠. For the three-variable Fenton–Karma model, we consider parameter sets 3, 4, 5, and 6 given by Fenton et al.²⁴ and reported in Table II. For the biophysically detailed ionic model of Ten tusscher et al. and Grandi et al., our simulations used the C++ implementations of the models provided by the authors, using $C_{m} = 1 μ F / {cm}^{2}$ for the membrane capacitance in both models.

For all cases, the domain is the interval $Ω = [0, 5]$ cm. The simulations use a mesh size $h = 31.25 μ$ m and a time step size $Δ t = 0.003653$ ms. An initial stimulus is applied at the center of the domain, $x \in [2.45, 2.55]$ cm, during the time interval $t \in [0.03, 1.03]$ ms. The conduction velocities are measured, as in the previous test case, by dividing the distance between two points x₁ and x₂ by the time interval between the activation times at these two selected locations. The computed conduction velocities of different ionic models are shown in Fig. 4 for relaxation times in the range $[0, 1]$ ms. Contrary to the results for the piecewise-linear ionic model, in this case, small and moderate values of the relaxation time act to enhance propagation.

FIG. 4.

View large Download slide

Dependence of the conduction velocity on the relaxation time for different ionic models: (left) Aliev–Panfilov; (center) Fenton–Karma; (right) ventricular Ten tusscher ’06 model and atrial Grandi ’11 model. An external stimulus is applied at the center of a one-dimensional domain and the activation times are recorded. For all models we used a time step size $Δ t = 0.003125$ ms and a mesh size $h = 31.25 μ$ m over a domain of length 5 cm. Unlike the linear case, the nonlinearities of the time derivative of the ionic currents increase the wavefront speed for small and moderate values of the relaxation time τ.

We also compare the computational time required by the parabolic and hyperbolic monodomain models in two and three spatial dimensions. Table III shows the number of iterations taken by the linear solvers (conjugate gradient preconditioned by successive over-relaxation) along with the relative computational times of the hyperbolic monodomain model (with respect to the parabolic model) for the reaction step and the diffusion step. In the reaction step, we solve the ionic model, and we also evaluate the ionic currents and their time derivatives. For this reason, we expect this step to take longer in the hyperbolic monodomain model. In fact, with simple ionic models, the amount of work required to evaluate the ionic currents and their time derivatives is almost doubled. This is natural because, for such simple models, the time derivative of the ionic currents do not use any quantities already computed in the evaluation of the ionic currents. By contrast, for biophysically detailed models, most of the computation is needed for the evaluation of the ionic currents, and many quantities can be reused for the evaluation of their time derivatives. This is reflected in Table III by a small increase (less than 7%) in the computational time of the hyperbolic model for the Ten Tusscher ’06 model. Unexpectedly, the solution of the linear systems take fewer iterations in the hyperbolic monodomain model. This is reflected by a speedup of about 20% of the computational time in the diffusion step, where we assembled the right hand side, solve the linear system, and update the variables Q and V.

TABLE III.

Comparison of the compuational cost for the parabolic (τ = 0 ms) and hyperbolic (⁠ $τ = 0.4$ ms) monodomain models for two- and three-dimensional problems in a serial computation. The negative values indicate relative slowdown of the hyperbolic model compared to the parabolic model, whereas the positive values represent relative speedup. In the reaction step, we measure the time needed to solve the ionic model and to evaluate the ionic currents and their time derivatives nodewise for 1000 time steps. In the diffusion step, we measure the time needed to form the right hand side, solve the linear system, and update the solutions, for 1000 time steps. The hyperbolic monodomain model can be solved using fewer linear solver iterations. This is reflected by the fact that the diffusion step is about 20% faster than in the parabolic monodomain model. On the other hand, the solution of the hyperbolic monodomain model requires the additional evaluation of the derivative of the ionic currents. For this reason, the reaction step is actually slower than in the parabolic monodomain model. Notice, however, the difference in the computational cost between the hyperbolic and parabolic reaction steps is smaller for the more complex ionic model. AP: Aliev–Panfilov ionic model; FK: Fenton–Karma ionic model; TP06: TenTusscher et al. ’06 ionic model.

2D	Max CG iterations		Reaction step cost	Diffusion step cost	Overall
100 × 100	τ = 0 ms	τ = 0.4 ms	$τ = 0.4 ms / τ = 0$ ms	$τ = 0.4 ms / τ = 0$ ms	$τ = 0.4 ms / τ = 0$ ms
AP	8	3	−54%	23%	−3%
FK	22	15	−62%	22%	−2%
TP06	11	5	−7%	15%	−3%

2D	Max CG iterations		Reaction step cost	Diffusion step cost	Overall
100 × 100	τ = 0 ms	τ = 0.4 ms	$τ = 0.4 ms / τ = 0$ ms	$τ = 0.4 ms / τ = 0$ ms	$τ = 0.4 ms / τ = 0$ ms
AP	8	3	−54%	23%	−3%
FK	22	15	−62%	22%	−2%
TP06	11	5	−7%	15%	−3%

3D	Max CG iterations		Reaction step cost	Diffusion step cost	Overall
10 × 10×10	τ = 0 ms	τ = 0.4 ms	$τ = 0.4 ms / τ = 0$	$τ = 0.4 ms / τ = 0$ ms	$τ = 0.4 ms / τ = 0$ ms
AP	43	16	−73%	48%	30%
FK	43	16	−49%	47%	25%
TP06	43	16	−6%	46%	15%

3D	Max CG iterations		Reaction step cost	Diffusion step cost	Overall
10 × 10×10	τ = 0 ms	τ = 0.4 ms	$τ = 0.4 ms / τ = 0$	$τ = 0.4 ms / τ = 0$ ms	$τ = 0.4 ms / τ = 0$ ms
AP	43	16	−73%	48%	30%
FK	43	16	−49%	47%	25%
TP06	43	16	−6%	46%	15%

C. Conduction velocity anisotropy

This section analyzes whether inductances affect the anisotropy in conduction velocity. We define the anisotropy ratio as the ratio between the conduction velocities in the transverse and longitudinal directions, and to obtain a simple estimate for this ratio, we consider plane-wave solutions propagating in these directions. We use the Fenton–Karma model with parameter set 3 and consider four values of the conductivity: $σ_{1} = 0.1$ mS/mm, $σ_{2} = 0.05$ mS/mm, $σ_{3} = 0.025$ mS/mm, and $σ_{4} = 0.0125$ mS/mm. Evaluating the wave speeds v₁, v₂, v₃, and v₄ corresponding to the conductivities σ₁, σ₂, σ₃, and σ₄, respectively, we show how the anisotropy ratio changes with respect to the relaxation time for several conductivity ratios. In fact, as shown in Table II, σ₁ and σ₄ are the typical longitudinal and transverse conductivities for this model. Considering the interval $Ω = [0, 5]$ cm, we use a mesh size $h = 25 μ$ m and a time step size $Δ t = 0.0025$ ms. Figure 5 (left) shows how the velocities at different σ compare to each other. To make the comparison more clear, we normalize the results with respect to the conduction velocities of the parabolic monodomain model, i.e., corresponding to τ = 0 ms. The results shown in Fig. 5 (left) should be understood as the percentage difference of the wave speed with respect to the parabolic case. It is clear that the curves are similar for the conductivities considered. On the other hand, the relaxation time at which the conduction velocity of the hyperbolic monodomain model matches the conduction velocity of the parabolic monodomain model decreases with smaller conductivities. Figure 5 (right) shows the relaxation time needed to maintain the same conduction velocity as in the parabolic monodomain model. In particular, for σ₄, we find that the relaxation time needed in the hyperbolic monodomain model to yield the same velocities as in the parabolic monodomain models is about $τ = 0.38$ ms. The differences in the velocities computed with $τ = 0.4$ ms and $τ = 0.38$ ms are smaller than 2%. This difference is typically smaller than the numerical error in the simulations. For this reason, in the following tests, our comparison will only use the value $τ = 0.4$ ms. Figure 5 (center) shows how the anisotropy ratio is influenced by the relaxation time. Although the ratio is not constant, the variation over the relaxation times considered never exceeds 5%. For the conductivity ratios 8:1 (σ₁:σ₄), 4:1 (σ₁:σ₃), and 2:1 (σ₁:σ₂), we find that the ratios between the conduction velocities in the transverse and in the longitudinal directions are approximately 1/ $\sqrt{8}$ (v₄/v₁), 1/2 (v₃/v₁), and 1/ $\sqrt{2}$ (v₂/v₁). This behavior is expected because the conduction velocity depends on the square root of the conductivity. We conclude that the effect of the inductances on the anisotropy ratio is negligible.

FIG. 5.

View large Download slide

Variation of the conduction velocity with respect to the conducitivity coefficients using the Fenton–Karma model and parameter set 3 in Table II. (Left) wave speed for different conductivities, normalized with respect to the parabolic conduction velocity. (Center) ratio between the conduction velocities as a function of the relaxation time. The variations are within 5% of the mean. (Right) relaxation time at which the hyperbolic conduction velocity matches the parabolic conduction velocity. Although the relaxation times at which we find the same wave speed for the parabolic and hyperbolic monodomain models are different for different conductivities, the error in the conduction velocity we commit considering the fixed value $τ = 0.4$ ms is smaller than 2% and therefore typically smaller than the spatial error.

D. Discretization error: A simple two-dimensional benchmark

It is important to be aware of the limitations of the numerical method used in the simulations. For example, spiral break up can occur easily if the wavefront is not accurately captured. In fact, numerical error can introduce spatial inhomogeneities that lead both to spiral wave formation and also to spiral break up. This numerical artifact disappears under grid refinement, and on sufficiently fine computational meshes, spiral wave break up results only from physical effects, such as conduction block.

To examine the role of spatial discretization on the system dynamics, we consider a square slab of tissue, $Ω = [0, 12] \times [0, 12]$ cm, and we use the Fenton–Karma model with the parameter set 3. Reentry is induced by an S1-S2 protocol. First, an external stimulus of unitary magnitude is applied in the left bottom corner, $Ω_{stim} = {x \in Ω : | | x | |_{1} \leq 1}$ at t = 0 ms. A second stimulus with the same amplitude and in the same region is applied after 300 ms. We consider two cases: (1) the fiber field is aligned with some of the mesh edges; and (2) the fiber field is not aligned with any edge in the mesh. Figures 6 (top) and 7 (top) demonstrate that this can be easily achieved by fixing the fiber field and rotating the mesh. Those figures show the fiber field in green on top of the computational grid. The red region at the bottom left of the domain is the stimulus region, $Ω_{stim} .$ In the first test (Fig. 6), the fibers are set to be orthogonal to the propagating front. Using 512 elements per side, which is equivalent to a mesh size h of approximately 234 μm (163 μm in the direction of wave propagation), the solutions obtained on the two meshes differ substantially. Large differences can also be found when we use 1024 elements per side, which is equivalent to a mesh size h of approximately 117 μm (82.5 μm in the direction of wave propagation). It is clear from Fig. 6 that whenever the fiber field is not aligned with the mesh, the conduction velocity is largely overestimated. The second test is similar to the first one, but the fibers are now rotated by 90°; see Fig. 7. It is clear that the error on the conduction velocity in the fiber direction is much smaller than the error in the transverse direction. In fact, the solutions for $h \approx 234 μ m$ and $h \approx 117 μ m$ show the greatest differences for transverse propagation.

FIG. 6.

View large Download slide

Effect of fiber direction and mesh orientation on propagation. (Top) fibers field (green), mesh orientation, and stimulus region (red); (Left) monodomain model; (Right) hyperbolic monodomain model. The domain size is 12 × 12 cm and the solution is evaluated at time t = 450 ms.

FIG. 7.

View large Download slide

Effect of fiber direction and mesh orientation on propagation. (Top) fibers field (green), mesh orientation, and stimulus region (red); (Left) monodomain model; (Right) hyperbolic monodomain model. The domain size is 12 × 12 cm and the solution is evaluated at time t = 360 ms.

Although similar results have already been reported in prior work,^49,61,66 it is nonetheless generally believed that a mesh size of the order of 200 μm is usually sufficient for cardiac computational electrophysiology.^23,48 This belief is supported by numerical simulations using a mesh size of 200 μm that show that the error on the wave speed in the fiber direction is less than 5%. On the other hand, the most interesting dynamics of the propagating front often will occur perpendicular to the fiber direction. For example, during the normal electrical activation of the ventricles, the signal spreads from the endocardium to the epicardium traveling across the ventricular wall, perpendicular to the fibers. The reduced conductivity in the transversal direction, usually taken to be about 8 times smaller, requires a finer resolution of the grid. As noted by Quarteroni et al.,⁶⁸ the required grid resolution to capture the transverse conduction velocity with an error smaller than 5% is about 25 μm. These results strongly indicate that to capture correctly the wave speeds, the mesh discretization needs to be about the size of a single cardiomyocyte, as also discussed by Hand and Griffith.³²

E. Effect of the relaxation time on spiral break up

Our next tests explore the effect of relaxation time on spiral wave break up. We consider a square slab of tissue, $Ω = [0, 12] \times [0, 12]$ cm. An initial stimulus is applied in the region y < 0.5 cm for 1 ms at t = 0 ms to generate a wave propagating in the y-axis. A second stimulus is then applied in the region ${x < 6 cm \land y < 7 cm}$ for 1 ms at time t = 320 ms to initiate a spiral wave. Spiral break up is easily obtained using the parameter set 3 in Table II, because of the steep action potential duration (APD) restitution curve. In this case, the back of the wave forms scallops, and when the turning spiral tries to invade these regions, it encounters refractory tissue and breaks.

Figure 8 shows the formation of the spiral wave in the case of isotropy. With $τ = 0.4$ ms and the Fenton–Karma model with parameter set 3, the conduction velocity of the hyperbolic monodomain model is very close to the conduction velocity of the parabolic monodomain model. In fact, the evolution of the spiral waves in the two cases is very similar.

FIG. 8.

View large Download slide

Comparison of the evolution of a spiral wave for the isotropic monodomain model (left) and the isotropic hyperbolic monodomain model (right) using the Fenton–Karma model with parameter set 3. The wavefronts (red) and the tails (light blue) represented by the variable Q are shown next to the transmembrane potential V. The evolution of the spiral wave is similar in the two cases. The domain is $[0, 12] \times [0, 12]$ cm, the mesh size is $h \approx 100 μ$ m, and the time step size $Δ t = 0.125$ ms. The fibers are aligned with the y-axis.

Figure 9 compares spiral break up obtained using the monodomain model and the hyperbolic monodomain model in the anisotropic case, with fibers aligned with the y-axis. The relaxation time for the hyperbolic monodomain model is chosen to be $τ = 0.4$ ms. As explained in Sec. IV C, although the conduction velocities in the transverse direction are different in the hyperbolic and parabolic models for $τ = 0.4$ ms, the difference is smaller than 2%. This difference is smaller than the spatial error and, for this reason, we consider the two cases to yield essentially the same anisotropy ratio. It is clearly shown in Fig. 9 that in the initial phase (up to t = 600 ms), the spirals are almost identical. After the first rotation, however, the spiral wave starts to break, and the relaxation time of the hyperbolic monodomain model shows quantitative differences in the form of the break up. Figure 9 shows the dynamics of $Q = \partial_{t} V$ ⁠. This variable can be used to define the fronts (red) and the tails (light blue) of the waves.

FIG. 9.

View large Download slide

Comparison of the evolution of a spiral break up for the monodomain model (left) and the hyperbolic monodomain model (right) using the Fenton–Karma model with parameter set 3. The wavefronts (red) and the tails (light blue) represented by the variable Q are shown next to the transmembrane potential V. Although the initiation sequence is the same and the spiral wave at t = 600 ms are similar, the evolution of the break up is different. The domain is $[0, 12] \times [0, 12]$ cm, the mesh size is $h \approx 200 μ$ m in the fiber direction, and $h \approx 50 μ$ m in the transverse direction, and the time step is $Δ t = 0.125$ ms. The fibers are aligned with the y-axis.

F. Atrial fibrillation

As a more complex application of the model, we use the hyperbolic monodomain model to simulate atrial fibrillation. Although a rigorous study of atrial fibrillation requires the use of a bidomain model, similar to the one proposed in Eqs. (11) and (12), the results given by the monodomain model can represent a reasonable approximation to the dynamics of the full bidomain model. For instance, Potse et al.⁶⁷ provide a comparison of the dynamics of the bidomain and monodomain models.

The anatomical geometry used in these simulations was based on the human heart model constructed by Segars et al.⁷⁶ using a 4D extended cardiac-torso (XCAT) phantom. From their data, we extracted and reconstructed a geometrical representation of the left atrium using SOLIDWORKS. We then used Trelis to generate a simplex mesh of the left atrium consisting of approximately 3.5 × 10⁶ elements. Our simulations use the Fenton–Karma model with parameter set 3.

A major challenge in modeling the atria is the definition of the fiber field. In work spanning the past 100 years, several studies have tried to characterize the fiber architecture of the atria,^62,63,73 but the structure of the muscle in the atria is so complex (as can be seen from the diagrams by Thomas⁸²) that developing a set of mathematical rules to reproduce a realistic fiber field is challenging. Our strategy for modeling the fiber field is to use the gradient of a harmonic function that satisfies prescribed boundary conditions. We divided the left atrium into several regions, and solved three Poisson problems: one for the left atrial appendage; one for the Bachmann’s bundle; and one for the remaining parts of the left atrium. In particular, for the left atrial appendage, we used a method similar to the one proposed by Rossi et al.⁶⁹ In the other regions, similar to work by Patelli et al.⁶⁴ and Krueger et al.,⁵⁰ we used the normalized gradient of the solution of the Poisson problem to reconstruct a fiber field similar to the one shown in the studies by Ho et al.^35,36 The approximate fiber field of the left atrium reconstructed using this strategy is shown in Fig. 10, where the colors represent the magnitude of the x-component of the fibers and are used only to highlight changes in the direction.

FIG. 10.

View large Download slide

(Left) Reconstructed fiber field of left atria based on the data in. (Right) activation times for the monodomain and for the hyperbolic monodomain model with $τ = 0.4$ ms using the Fenton–Karma model with parameter set 3.

To initiate atrial fibrillation, we use an S1-S2 stimulation protocol that is applied at the junction with the inter-atrial band, as described by Colman.¹⁹ The S1 stimulus is applied with a cycle length of 350 ms followed by a short S2 stimulus with a cycle length of 160 ms. Figure 10 (right) shows that the initial electrical activation times are similar for the monodomain and hyperbolic monodomain models. Because this set of parameters gives a velocity of about 45 cm/s, roughly 35% slower than the expected conduction velocity,³⁴ it takes longer for the atrium to be fully activated. When the second stimulus is applied, spiral waves form in the hyperbolic monodomain system but not in the parabolic monodomain model; see Fig. 11 at 420 ms. Spiral waves are generated for τ = 0 ms only after the third stimulus is applied. Figure 11 also shows the fronts (red) and tails (light blue) of the waves, as indicated by Q.

FIG. 11.

View large Download slide

Comparison of spiral break up in the left atria using a fast pacing S1-S2 stimulation protocol using the Fenton–Karma model with parameter set 3. On the left we show the solution of the monodomain model, relaxation time τ = 0 ms, and on the right we show the solution of the hyperbolic monodomain model, relaxation time $τ = 0.4$ ms. Using the hyperbolic monodomain model spiral waves form after second stimulus. The parabolic monodomain model requires 3 stimuli before spiral waves are formed. The top row show the solution for the transmembrane potential, and the bottom row shows the wavefronts (red) and the tails (light blue) represented by the variable Q.

G. Virtual electrodes

The reduction of the bidomain model to the monodomain model is not valid if the extracellular and intracellular compartments have different anisotropy ratios. As shown in experimental studies,^45,85 when a current is supplied by an extracellular electrode, adjacent areas of depolarization and hyperpolarization are formed in the case of unequal anisotropy ratios. Applying a stimulation current $I_{e, stim}$ in the extracellular compartment, the region of hyperpolarization near a cathode (⁠ $I_{e, stim} < 0$ ⁠) is called a virtual anode, while the region of depolarization near an anode (⁠ $I_{e, stim} > 0$ ⁠) is called a virtal cathode. The existence of virtual cathodes and anodes is important to understand the four mechanisms of cardiac stimulation and to study the mechanisms of defribillation. In this section, we focus only on a unipolar cathodal stimulation to investigate the cathode formation mechanism in the hyperbolic bidomain model. We refer the reader interested in this virtual electrode phenomenon to the more detailed numerical studies of Wikswo and Roth⁸⁶ and Colli Franzone et al.¹⁸

Following Sepulveda et al.,⁷⁷ we consider the domain $Ω = [- 2, 2] cm \times [- 0.8, 0.8] cm$ ⁠, and consider this to correspond to the epicardial surface. We apply a unipolar cathodal current in the region $Ω_{stim} = [- 0.5, 0.5] mm \times [- 0.1, 0.1] mm$ ⁠. We use the Fenton–Karma model with parameter set 3, and the corresponding nondimensionalized current stimulus amplitude is –100. We choose $σ_{e}^{f} = 1.5448$ ⁠, $σ_{e}^{s} = σ_{e}^{n} = 1.0438$ ⁠, and $σ_{i}^{f} = 2.3172$ ⁠, $σ_{i}^{s} = σ_{i}^{n} = 0.2435$ ⁠, with the fibers aligned with the x-axis.

Figure 12 shows the pattern formed by the transmembrane potential 2 ms after the extracellular stimulus is applied. The red regions are areas of depolarization and the blue regions are the areas of hyperpolarization. The depolarization and hyperpolarization regions are affected by the choice of the extracellular and intracellular relaxation times. For this test, we consider the relaxation times $τ_{i}$ and $τ_{e}$ to take the values 0, 0.2, and 0.4 ms. Note that if $τ_{i} = τ_{e}$ ⁠, the hyperbolic bidomain equations (13) and (14) reduce to the hyperbolic-elliptic system of Eqs. (15) and (16). Although this test shows the ability of the hyperbolic bidomain model to reproduce the virtual electrode phenomenon, the influence of the relaxation times in the stimulation remains unclear and necessitates further investigation.

FIG. 12.

View large Download slide

Transmembrane potential patterns formed by the hyperbolic bidomain model with the Fenton–Karma model at 2 ms after a unipolar cathode stimulus is supplied in the extracellular compartment. The red regions are the areas of depolarization and blue regions are the areas of hyperpolarization. The hyperbolic bidomain model shows a behavior qualitatively similar to the conventional bidomain model shown in the study by Ref. 86, but the details are influenced by the extracellular and intracellular relaxation times.

V. CONCLUSIONS

Local perturbation in the standard bidomain and monodomain models propagates with an infinite speed. To correct this unrealistic feature, we developed a hyperbolic version of the bidomain model, and in the case that intracellular and extracellular conductivity tensors have equal anisotropy ratios, we reduced this model to a hyperbolic monodomain model. Our derivation relies on a Cattaneo-type model for the fluxes, described by Eqs. (1) and (2). The Cattaneo-type fluxes are equivalent to introducing self-inductance effects, as shown by the schematic diagram in Fig. 1. As shown in Appendix A, relaxation times introduced in the Cattaneo fluxes represent the ratio between the inductances and the resistances of the extracellular and intracellular compartments. Although the hyperbolic monodomain model reduces to the classical parabolic model in the case that the relaxation times of both the intracellular and extracellular compartments are zero, the models differ if at least one of these relaxation times is nonzero.

Although hyperbolic bidomain and monodomain models do not appear to have been previously discussed in the literature, the work of Hurtado et al.³⁹ does address the problem of propagation with infinite speed. Their approach uses a nonlinear model for the fluxes based on porous medium assumptions. Because of the nonlinear nature of the fluxes, the porous medium approach to cardiac electrophysiology is difficult to analyze, even for simple linear reactions. Hurtado et al.³⁹ use the simplified ionic model proposed by Aliev and Panfilov.³ It would be interesting to see the extension to biophysically detailed models. Moreover, it is clear that the porous medium approach is not incompatible with the theory developed here. In fact, the two models could be combined using a nonlinear Cattaneo-type model for the fluxes.

As discussed by King et al.,⁴³ abnormal conduction velocities may contribute to arrhythmogenesis. We have studied how the conduction velocities of the propagating action potential change with respect to the relaxation time of the Cattaneo-type fluxes. In the case of the piecewise-linear model of McKean,⁵⁵ it is possible to find an analytical expression for the conduction velocity of the hyperbolic monodomain model. (Unfortunately, the wave speed of this bistable model has been reported incorrectly in some prior studies.^1,56) We have shown that for the hyperbolic monodomain model, signals cannot travel faster than the characteristic propagation speed of the medium. Our numerical results match the theoretical predictions for different values of the excitability parameter. In this simple case, conduction velocity is a monotone decreasing function of relaxation time. Additionally, we have shown that, for this linear case, discretization errors are lower in the hyperbolic model.

Using simple nonlinear ionic models,^3,25 however, we found that the relationship between the relaxation time and the conduction velocity is not always monotone. In fact, for small and moderate values of the relaxation time, we found that waves propagate faster than in the standard monodomain model. With these ionic models, a maximum conduction velocity is found at an optimal relaxation time. At small relaxation times, the second derivative terms are relatively small. Consequently, this phenomenon likely results from the time derivative of the ionic currents. For larger values of the relaxation times, inertial effects dominate, and the conduction velocity monotonically decreases to zero. This implies that we can find a value of the relaxation time for which the action potential propagates at the same velocity as in the standard monodomain model. This fact allowed us to directly investigate the influence of inductances in the monodomain model. In the biophysically detailed ionic models for ventricular myocytes⁸¹ and atrial myocytes³⁰ considered here, a similar behavior can be found, although only for very small values of the relaxation time. For this reason, we compared the parabolic and hyperbolic monodomain models in two and three spatial dimensions using the Fenton–Karma model.²⁵

Changing the type of equation may have important consequences in terms of the stability of the numerical approximations. We did not encounter any substantial numerical difficulties beyond those already present in the standard parabolic models. On the contrary, we have shown that the linear system is easier to solve and the solutions are more accurate when the hyperbolic system is used. Moreover, we believe that the solution of the parabolic monodomain and bidomain model already have many of the numerical difficulties usually associated with hyperbolic systems: the propagating action potential in the parabolic models already have sharp fronts which represent one of the main difficulties in hyperbolic systems. Additionally, the hyperbolic monodomain and bidomain models are damped wave equations, where the damping term is dominant. In conclusion, the introduction of hyperbolic terms to the monodomain and bidomain models does not appear to add substantial numerical difficulties.

We performed several qualitative comparisons between the parabolic and the hyperbolic monodomain models. In particular, we considered the case of spiral break up resulting from the conduction block caused by a steep APD restitution curve. Using parameter set 3 of the Fenton–Karma model, the back of the wave easily forms a series of indentations or scallops. When a spiral wave invades this after turning, it encounters refractory tissue that causes the conduction block, leading to a breakup of the wave. This behavior, already observed for the standard monodomain model and explained by Fenton et al.,²⁴ is still present in the hyperbolic monodomain model. We showed spiral break up in a simple two-dimensional test case and also in an anatomically realistic three-dimensional model of the left atrium. We reconstructed on the left atria a fiber field qualitatively in accordance with the data reported by Ho et al.,³⁵ Krueger et al.,⁵⁰ and Pashakhanloo et al.⁶³ We applied a fast-pacing S1-S2 stimulation protocol to the left atrium to induce spiral waves and spiral break up. Although the initial activation sequences were similar, for the hyperbolic monodomain model, spiral waves appeared already after the second stimulus was applied, whereas for the parabolic monodomain at least three stimuli were needed.

We also used a simple two-dimensional benchmark to investigate how the alignment of the fiber direction with respect to the elements of the finite element mesh influences wave propagation. Under mesh refinement, spatial discretization errors become smaller and smaller, but the common assumption that a grid resolution of 200 μm is sufficient to correctly capture the conduction velocity^48,49 seems to be overly optimistic. In our tests, we analyzed the propagation of a simple wave using an anisotropic conductivity tensor. Fixing the fiber field and rotating the mesh, we were able to show that using a mesh size of about 234 μm (163 μm in the direction of wave propagation), the propagation was greatly influenced by the mesh orientation relative to the fiber alignment. Even on a finer mesh with a mesh size of about 117 μm (82.5 μm in the direction of wave propagation), visible differences between the propagation patterns obtained with the different orientations remained. These tests highlight the fact that if the transverse conductivity coefficients are taken to be about 8 times smaller than the longitudinal conductivity, as is typically done in practice,^9,16,38,70 the mesh size needs to be much smaller than 200 μm in order to yield resolved dynamics. Our results are in accordance with the convergence study by Quarteroni et al.⁶⁸ for propagation in the transverse direction, in which the authors showed that the mesh size in the transverse direction needs to be about 25 μm. Because such mesh resolutions require elements only slightly larger than the actual cardiac cells, it is natural to ask whether the use of a continuum model is even appropriate: for a similar computational cost, it could be possible to construct a discrete model of the heart. Although the difficulties in representing correctly the conduction velocities are well documented in the literature^49,65 and have been thoroughly analyzed by Pezzuto et al.,⁶⁶ the focus generally seems to have been on the conduction in the fiber direction, as is clear from the choice of the three-dimensional benchmark test proposed by Niederer et al.⁶¹ Other solution strategies for the monodomain and bidomain models using adaptive mesh refinement^17,47 and high-order elements⁴ have been proposed, but their application so far appears to be relatively limited.

Finally, we showed that the hyperbolic bidomain model can capture the virtual electrode phenomenon observed in the standard model. By stimulating the epicardial surface using a unipolar cathodal current, we have shown how the pattern of the transmembrane potential is influenced by the intracellular and extracellular relaxation times. Nonetheless, a more detailed study on how the relaxation times affect the virtual electrodes is needed to fully understand the hyperbolic bidomain model. Comparing the stimulation patterns to experimental data could even serve to calibrate the magnitude of the relaxation times in the hyperbolic bidomain model.

Despite significant progress in models of cardiac electrophysiology, the complex effects of nonlinearity and heterogeneity in the heart are far from being fully understood. We have shown that the nonlinearities of the underlying physics can give unexpected results, in contradiction to the linear case. Inductances in tissue propagation could play an important role in cardiac electrical dynamics, especially close to the wavefronts, where the fast currents responsible for the initiation of the action potential can give a small but non-negligible contribution. The main phenomenon we observed in the hyperbolic model is the enhancement of the conduction velocity at small relaxation times because of the presence of the time derivative of the ionic currents. This phenomenon would allow electric signals to propagate at the same speed with lower conductance as compared to the standard models. Future experimental studies are needed to confirm the importance of these effects.

ACKNOWLEDGMENTS

The authors gratefully acknowledge research support from NIH Award No. HL117063 and from NSF Award No. ACI 1450327. The human atrial model used in this work was kindly provided by Professor W. P. Segars. We also thank Professor C. S. Henriquez and Professor J. P. Hummel for extended discussions on atrial modeling. We also gratefully acknowledge support by the libMesh developers in aiding the development of the code used in the simulations reported herein.

APPENDIX A: CABLE EQUATION WITH INDUCTANCES

We follow the derivation of the cable equation by Keener and Sneyd.⁴² As shown in Fig. 1, we assume there exist inductance effects in both the extra- and intracellular axial currents. We have

V_{e} (x) - V_{e} (x + Δ x) = - L_{e} \partial_{t} I_{e} Δ x - R_{e} I_{e} Δ x,

(A1)

V_{i} (x) - V_{i} (x + Δ x) = - L_{i} \partial_{t} I_{i} Δ x - R_{i} I_{i} Δ x,

(A2)

where $I_{e}$ and $I_{i}$ are the extracellular and intracellular axial currents, respectively. Dividing by $Δ x$ and taking the limit as $Δ x \to 0$ ⁠, we find

\partial_{x} V_{e} = - L_{e} \partial_{t} I_{e} - R_{e} I_{e},

(A3)

\partial_{x} V_{i} = - L_{i} \partial_{t} I_{i} - R_{i} I_{i} .

(A4)

The minus sign on the right hand sides is a convention that ensures that positive charges flow from left to right. Applying Kirchhoff’s current law, we have

I_{e} (x) - I_{e} (x + Δ x) + I_{t} Δ x = 0,

(A5)

I_{i} (x) - I_{e} (x + Δ x) - I_{t} Δ x = 0,

(A6)

which, in the limit $Δ x \to 0$ ⁠, becomes

I_{t} = \partial_{x} I_{e} = - \partial_{x} I_{i} .

(A7)

Differentiating Eqs. (A3) and (A4) with respect to x, and Eq. (A7) with respect to t, we find that for the extracellular compartment

\partial_{x x} V_{e} = - L_{e} \partial_{x t} I_{e} - R_{e} \partial_{x} I_{e},

(A8)

\partial_{x t} I_{e} = \partial_{t} I_{t},

(A9)

and for the intracellular compartment,

\partial_{x x} V_{i} = - L_{i} \partial_{x t} I_{i} - R_{i} \partial_{x} I_{i},

(A10)

\partial_{x t} I_{e} = - \partial_{t} I_{t} .

(A11)

Substituting the mixed derivative of the extracellular and intracellular currents, we find

\partial_{x x} V_{e} = - L_{e} \partial_{t} I_{t} - R_{e} I_{t}, \partial_{x x} V_{i} = L_{i} \partial_{t} I_{t} + R_{i} I_{t} .

(A12)

Defining $V = V_{i} - V_{e}$ ⁠, and taking the difference between the intracellular and the extracellular equations, we can write a single equation for V

\partial_{x x} V = (L_{i} + L_{e}) \partial_{t} I_{t} + (R_{i} + R_{e}) I_{t} .

(A13)

Recalling that

I_{t} = χ (C_{m} \partial_{t} V + I_{ion}),

(A14)

we finally find the hyperbolic form of the cable equation

τ C_{m} \frac{\partial^{2} V}{\partial t^{2}} + C_{m} \frac{\partial V}{\partial t} - \frac{\partial}{\partial x} (D \frac{\partial V}{\partial x}) = - I_{ion} - τ \frac{\partial I_{ion}}{\partial t},

where we define the conductivity as $D = 1 / χ (R_{i} + R_{e})$ and the relaxation time by

τ = \frac{L_{i} + L_{e}}{R_{i} + R_{e}} .

(A15)

APPENDIX B: EXACT SOLUTION FOR THE PIECEWISE-LINEAR BISTABLE MODEL

Consider the simplified piecewise-linear bistable model

\begin{matrix} I_{ion} (V) = k V - k [V_{2} H (V - V_{1}) + V_{0} (1 - H (V - V_{1}))] . \end{matrix}

(B1)

The nondimensional form of Eq. (B1) is the simplified model proposed by McKean,⁵⁵

{\hat{I}}_{ion} (\hat{V}) = \hat{V} - H (\hat{V} - α) .

(B2)

Using (B1) in the hyperbolic monodomain model (21) and considering only one spatial dimension, we have

{\begin{matrix} μ \frac{\partial^{2} \hat{V}}{\partial {\hat{t}}^{2}} + (1 + μ) \frac{\partial \hat{V}}{\partial \hat{t}} - \frac{\partial^{2} \hat{V}}{\partial {\hat{x}}^{2}} + \hat{V} = 0, & \hat{V} < α, \\ μ \frac{\partial^{2} \hat{V}}{\partial {\hat{t}}^{2}} + (1 + μ) \frac{\partial \hat{V}}{\partial \hat{t}} - \frac{\partial^{2} \hat{V}}{\partial {\hat{x}}^{2}} + \hat{V} = 1, & \hat{V} > α, \end{matrix}

(B3)

where we have introduced the nondimensional variables

\hat{t} = \frac{1}{T} \hat{t}, \hat{x} = \frac{1}{L} x, \hat{V} = \frac{V - V_{0}}{V_{2} - V_{0}} .

In particular, we define

T = \frac{C_{m}}{k}, L = \sqrt{\frac{σ}{k χ}}, μ = \frac{τ k}{C_{m}},

where μ is a nondimensional number that characterizes the effect of the inductances in the system. Specifically, μ is the ratio between the relaxation time of the system and the characteristic time of the reactions: the larger μ, the more important the effects of the inductances become. Introducing the change of coordinates $z = x - c t,$ such that $U (z) = U (x - c t) = \hat{V} (x, t)$ ⁠, the system (B3) is transformed into

{\begin{matrix} (c^{2} μ - 1) U_{z z} - c (1 + μ) U_{z} + U = 0, & z > 0, \\ (c^{2} μ - 1) U_{z z} - c (1 + μ) U_{z} + U = 1, & z < 0. \end{matrix}

(B4)

Consider first the case z > 0. Using $γ = c^{2} μ - 1$ and $β = - c (1 + μ)$ ⁠, this equation reads

γ V ″ + β V' + V = 0

for which the roots of the characteristic polynomial are

λ_{\pm} = \frac{- β \pm \sqrt{β^{2} - 4 γ}}{2 γ},

(B5)

so that the solution is $U (z) = A_{+} e^{λ_{+} z} + A_{-} e^{λ_{-} z}$ ⁠. It is easy to verify that the case z < 0 has the solution $U (z) = B_{+} e^{λ_{+} z} + B_{-} e^{λ_{-} z} + 1$ ⁠. The global solution, therefore, is

U (z) = {\begin{matrix} A_{+} e^{λ_{+} z} + A_{-} e^{λ_{-} z}, & z > 0, \\ B_{+} e^{λ_{+} z} + B_{-} e^{λ_{-} z} + 1, & z < 0. \end{matrix}

(B6)

For (B6) to represent a propagating front, it is necessary that the solution is real and bounded for any value of z. This implies that either $λ_{-} < 0 < λ_{+}$ or $λ_{+} < 0 < λ_{-}$ ⁠. Requiring $β^{2} - 4 γ > 0$ ⁠, we find the condition $γ < β^{2} / 4$ ⁠, which is always satisfied for the model (B4) because $β^{2} - 4 γ = c^{2} {(1 - μ)}^{2} + 4 > 0.$ If c > 0 then $β < 0$ ⁠. Then $λ_{-} > 0$ and $λ_{+} < 0$ if

c < \sqrt{\frac{1}{μ}} .

(B7)

Imposing $U (\infty) = 0$ ⁠, we find $A_{-} = 0$ ⁠, and imposing $U (- \infty) = 1$ ⁠, we find $B_{+} = 0$ ⁠, so that

U (z) = {\begin{matrix} A_{+} e^{λ_{+} z}, & z > 0, \\ B_{-} e^{λ_{-} z} + 1, & z \leq 0. \end{matrix}

(B8)

To ensure the continuity of the solution at z = 0, we fix $U (0) = α$ ⁠, so that $A_{+} = α = 1 + B_{-}$ and

U (z) = {\begin{matrix} α e^{λ_{+} z}, & z > 0, \\ (α - 1) e^{λ_{-} z} + 1, & z \leq 0. \end{matrix}

(B9)

The first derivative of U is

U' (z) = {\begin{matrix} α λ_{+} e^{λ_{+} z}, & z > 0, \\ (α - 1) λ_{-} e^{λ_{-} z}, & z \leq 0. \end{matrix}

(B10)

Imposing the continuity of $U'$ at z = 0, we find $α λ_{+} = (α - 1) λ_{-}$ and therefore

\frac{γ}{β^{2}} = \frac{(α^{2} - α)}{{(2 α - 1)}^{2}} .

(B11)

Using the definitions of γ and β in (B11), we finally obtain an expression for the speed c in terms of μ and α

c = \frac{(1 - 2 α)}{\sqrt{μ + (α - α^{2}) {(μ - 1)}^{2}}} \leq \sqrt{\frac{1}{μ}} .

(B12)

Notice that the characteristic propagation speed is bounded from above by $\sqrt{1 / μ}$ ⁠.

APPENDIX C: FINITE ELEMENT DISCRETIZATION

The common practice to solve the coupled nonlinear system of cardiac electrophysiology is to split the monodomain (or bidomain) model from the ionic model. Motivated by the work of Munteanu et al.,⁵⁸ we will also decouple the two systems. For brevity, we show the discretization only for the monodomain model with a first-order time integrator, because first-order time integrators are still widely used in the community. The same approach can be used to discretize the hyperbolic bidomain model. Given the subinterval $[t^{n}, t^{n + 1}]$ ⁠, with $t^{0} = 0$ ⁠, we define two separate subproblems, one for the ionic model and one for the monodomain equation connected by the initial conditions. In the first step, we solve the ionic model (25) for $w^{n + 1}$ ⁠, and we compute the ionic currents using the updated values of the state variables. The ionic currents computed in this way are then used to solve the monodomain system (22)–(23).

As shown in Spiteri and Dean,⁷⁹ the most efficient numerical algorithm to solve the stiff ordinary differential equations of the ionic model strongly depends on the model considered. We use the forward Euler method for the simplified ionic models of Aliev–Panfilov and Fenton–Karma. For the ventricular ionic model of Tusscher and Panfilov,⁸¹ we use the Rush–Larsen method.⁷¹ The atrial ionic model of Grandi et al.³⁰ has a stricter restriction on the time step size. Consequently, we use the Rush–Larsen method in combination to the backward Euler method for linear equations. The time step size is defined as $Δ t = t^{n + 1} - t^{n}$ ⁠. To obtain a second-order time discretization scheme, the ionic model can be solved using the explicit trapezoidal method (i.e., Heun’s method), which is a second-order accurate strong stability preserving Runge–Kutta method.

Once the solution of the ionic model has been found, we solve the monodomain model using a low-order finite element discretization. Because (22) is a simple ordinary differential equation, the computational costs for solving (21) or the system (22) and (23) is comparable. Denoting by ${(v, w)}_{Ω} = \int_{Ω} v w$ the $L^{2} (Ω)$ inner product and using the boundary conditions (24), the Galerkin approximation of Eqs. (22) and (23) is to find $V^{h}, Q^{h} \in S^{h}$ such that

{(\frac{\partial V^{h}}{\partial t}, ψ^{h})}_{Ω} = {(Q^{h}, ψ^{h})}_{Ω},

(C1)

\begin{matrix} {(τ C_{m} \frac{\partial Q^{h}}{\partial t} + C_{m} Q^{h}, ϕ^{h})}_{Ω} + {(D \nabla V, \nabla ϕ^{h})}_{Ω} \\ = {(I_{ion}^{h}, ϕ^{h})}_{Ω} + {(τ \frac{\partial I_{ion}^{h}}{\partial t}, ϕ^{h})}_{Ω} \end{matrix}

(C2)

for all $ψ^{h}, ϕ^{h} \in S^{h}$ ⁠, where

S^{h} = {v^{h} \in C^{0} (\bar{Ω}) : {v^{h} |}_{K} \in P^{1} (K), \forall K \in T^{h}} .

(C3)

In practice, we look for a continuous solution that is linear in every simplex element K in the triangulation $T^{h}$ of Ω. Introducing the basis of nodal shape functions ${N_{A}}_{A = 1}^{M},$ with $M = dim (S^{h})$ ⁠, the solution fields are discretized as

V^{h} = \sum_{A = 1}^{M} N_{A} V_{A}, Q^{h} = \sum_{A = 1}^{M} N_{A} Q_{A} .

Equations (C1) and (C2) yield the matrix system

\dot{V} = Q,

(C4)

C_{m} M (τ \dot{Q} + Q) + KV = F + τ L .

(C5)

We use the half-lumping scheme proposed by Pathmanathan et al.,⁶⁵ so that

F = MI, L = MJ

and

\dot{V} = Q,

(C6)

C_{m} M_{L} (τ \dot{Q} + Q) + KV = M (I + τ J),

(C7)

where

[M_{A B}] = {(N_{A}, N_{B})}_{Ω}, [K_{A B}] = {(\nabla N_{A}, D \nabla N_{B})}_{Ω},

$I$ and $J$ are the nodal evaluation of $I_{ion}$ and $\partial I_{ion} / \partial t$ ⁠, and $M_{L}$ is the lumped mass matrix obtained by row summation of the mass matrix. Using a simple first-order implicit-explicit time integrator, the fully discrete system of equations reads

V^{n + 1} = V^{n} + Δ t Q^{n + 1},

(C8)

\begin{matrix} C_{m} (τ + Δ t) M_{L} Q^{n + 1} + Δ t K V^{n + 1} = Δ t M (I^{*} + τ J^{*}) + τ C_{m} M_{L} Q^{n}, \end{matrix}

(C9)

where $I^{*}$ and $J^{*}$ indicate the values of the ionic currents, and of their derivatives, evaluated using the updated values for $w^{n + 1}$ obtained by solving the ionic model. Introducing (C8) in (C9), we arrive at a single system for $Q^{n + 1}$ ⁠,

\begin{matrix} [C_{m} (τ + Δ t) M_{L} + Δ t^{2} K] Q^{n + 1} \\ = τ C_{m} M_{L} Q^{n} - Δ t K V^{n} + Δ t M (I^{*} + τ J^{*}) . \end{matrix}

(C10)

The time derivative of the ionic currents, $J^{*}$ ⁠, is approximated nodally using the chain rule as

\begin{matrix} \frac{\partial I_{ion}}{\partial t} (t^{n}) \approx \frac{\partial I_{ion}}{\partial V} (V^{n}, w^{n + 1}) Q^{n} + \frac{\partial I_{ion}}{\partial w} (V^{n}, w^{n + 1}) \cdot g (V^{n}, w^{n + 1}), \end{matrix}

(C11)

where the function $g (V, w)$ is the right hand side of the ionic model system defined in Eq. (25). After solving (C10), we update the voltage Eq. (C8).

The scheme used to obtain (C8) and (C9) is the simplest method—ARS(1,1,1), from Ascher, Ruuth, and a series of implicit-explicit (IMEX) Runge–Kutta (RK) algorithms that are popular for hyperbolic systems.^5,11 For a second-order time integrator, we use the H-CN(2,2,2) scheme,¹⁰ which combines Heun’s method for the explicit part and the Crank–Nicholson method for the implicit part.

We have not experienced any difference in time step size restriction between the parabolic and hyperbolic models. In fact, we have used the IMEX-RK method where only the ionic currents and their time derivatives are treated explicitly in both parabolic and hyperbolic models. The time step size restriction is dictated only by the ionic model. Nonetheless, the time step size should be chosen to accurately capture the propagation of the action potential. Denoting with v the conduction velocity and with h_m the smallest element size, we chose our time step size such that the CFL condition $Δ t \leq h_{m} / v$ holds.

From a purely numerical perspective, in the first-order scheme, the dissipative nature of the backward Euler method can be used to remove numerical oscillations, while the second-order Crank–Nicholson method will keep such oscillations. A dissipative second-order time integrator could be used in place of the Crank–Nicholson method. In any case, the physical dissipation of the hyperbolic monodomain model reduces such artifacts for first- and second-order schemes.

One of the main difficulties arising from hyperbolic systems is their tendencies to form discontinuities for which numerical approximations typically develop spurious oscillations. Although the standard monodomain and bidomain models are parabolic systems, if the wavefront is not accurately resolved then the upstroke of the action potential can act as a discontinuity. For this reason, the numerical methods used for cardiac electrophysiology should be already suitable for the hyperbolic systems considered here.

References

1.

H. A.

Abdusalam

, “

Analytic and approximate solutions for Nagumo telegraph reaction diffusion equation

,”

Appl. Math. Comput.

157

(

2

),

515

–

522

(

2004

).

https://doi.org/10.1016/j.amc.2003.08.050

Google Scholar

Crossref

2.

J.

Ahrens

,

B.

Geveci

,

C.

Law

,

C.

Hansen

, and

C.

Johnson

,

36-paraview: An end-user tool for large-data visualization, Visualization Handb.

(

Butterworth-Heinemann

,

2015

), Chap. 36, pp.

717

–

731

.

https://doi.org/10.1016/B978-012387582-2/50038-1

Google Scholar

Crossref

3.

R. R.

Aliev

and

A. V.

Panfilov

, “

A simple two-variable model of cardiac excitation

,”

Chaos, Solitons Fractals

7

(

3

),

293

–

301

(

1996

).

https://doi.org/10.1016/0960-0779(95)00089-5

Google Scholar

Crossref

4.

C. J.

Arthurs

,

M. J.

Bishop

, and

D.

Kay

, “

Efficient simulation of cardiac electrical propagation using high order finite elements

,”

J. Comput. Phys.

231

(

10

),

3946

–

3962

(

2012

).

https://doi.org/10.1016/j.jcp.2012.01.037

Google Scholar

Crossref

PubMed

5.

U. M.

Ascher

,

S. J.

Ruuth

, and

R. J.

Spiteri

, “

Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations

,”

Appl. Numer. Math.

25

(

2–3

),

151

–

167

(

1997

).

https://doi.org/10.1016/S0168-9274(97)00056-1

Google Scholar

Crossref

6.

S.

Balay

,

S.

Abhyankar

,

M. F.

Adams

,

J.

Brown

,

P.

Brune

,

K.

Buschelman

,

L.

Dalcin

,

V.

Eijkhout

,

W. D.

Gropp

,

D.

Kaushik

,

M. G.

Knepley

,

L. C.

McInnes

,

K.

Rupp

,

B. F.

Smith

,

S.

Zampini

, and

H.

Zhang

, “

PETSc users manual

,”

Technical Report No. ANL-95/11-Revision 3.6

, Argonne National Laboratory (

2015a

).

Google Scholar

7.

S.

Balay

,

S.

Abhyankar

,

M. F.

Adams

,

J.

Brown

,

P.

Brune

,

K.

Buschelman

,

L.

Dalcin

,

V.

Eijkhout

,

W. D.

Gropp

,

D.

Kaushik

,

M. G.

Knepley

,

L. C.

McInnes

,

K.

Rupp

,

B. F.

Smith

,

S.

Zampini

, and

H.

Zhang

, http://www.mcs.anl.gov/petsc for “PETSc Web page” (

2015b

).

8.

S.

Balay

,

W. D.

Gropp

,

L. C.

McInnes

, and

B. F.

Smith

, “

Efficient management of parallelism in object oriented numerical software libraries

,” in

Modern Software Tools in Scientific Computing

, edited

E.

Arge

,

A. M.

Bruaset

, and

H. P.

Langtangen

(

Birkhäuser Press

,

1997

), pp.

163

–

202

.

Google Scholar

Crossref

9.

A. P.

Benson

,

O.

Bernus

,

H.

Dierckx

,

S. H.

Gilbert

,

J. P.

Greenwood

,

A. V.

Holden

,

K.

Mohee

,

S.

Plein

,

A.

Radjenovic

,

M. E.

Ries

et al., “

Construction and validation of anisotropic and orthotropic ventricular geometries for quantitative predictive cardiac electrophysiology

,” in

Interface Focus

(

The Royal Society, United Kingdom

,

2010

) p.

rsfs20100005

.

Google Scholar

Crossref

10.

S.

Boscarino

,

F.

Filbet

, and

G.

Russo

, “

High order semi-implicit schemes for time dependent partial differential equations

,”

J. Sci. Comput.

68

(

3

),

975

–

1001

(

2016

).

https://doi.org/10.1007/s10915-016-0168-y

Google Scholar

Crossref

11.

S.

Boscarino

and

G.

Russo

, “

On a class of uniformly accurate imex Runge–Kutta schemes and applications to hyperbolic systems with relaxation

,”

SIAM J. Sci. Comput.

31

(

3

),

1926

–

1945

(

2009

).

https://doi.org/10.1137/080713562

Google Scholar

Crossref

12.

A.

Bueno-Orovio

,

E. M.

Cherry

, and

F. H.

Fenton

, “

Minimal model for human ventricular action potentials in tissue

,”

J. Theor. Biol.

253

(

3

),

544

–

560

(

2008

).

https://doi.org/10.1016/j.jtbi.2008.03.029

Google Scholar

Crossref

PubMed

13.

A.

Bueno-Orovio

,

D.

Kay

,

V.

Grau

,

B.

Rodriguez

, and

K.

Burrage

, “

Fractional diffusion models of cardiac electrical propagation: Role of structural heterogeneity in dispersion of repolarization

,”

J. R. Soc. Interface

11

(

97

),

20140352

(

2014

).

https://doi.org/10.1098/rsif.2014.0352

Google Scholar

Crossref

PubMed

14.

C.

Cattaneo

, “

Sulla conduzione del calore

,”

Atti Semin. Mat. Fis. Univ. Modena

3

,

83

–

101

(

1948

).

Google Scholar

15.

D. E.

Clapham

and

L. J.

Defelice

, “

Membrane biolow small signal impedance of heart cell membranes

,”

J. Membr. Biol.

67

,

63

–

71

(

1982

).

https://doi.org/10.1007/BF01868648

Google Scholar

Crossref

PubMed

16.

R. H.

Clayton

,

O.

Bernus

,

E. M.

Cherry

,

H.

Dierckx

,

F. H.

Fenton

,

L.

Mirabella

,

A. V.

Panfilov

,

F. B.

Sachse

,

G.

Seemann

, and

H.

Zhang

, “

Models of cardiac tissue electrophysiology: Progress, challenges and open questions

,”

Prog. Biophys. Mol. Biol.

104

(

1

),

22

–

48

(

2011

).

https://doi.org/10.1016/j.pbiomolbio.2010.05.008

Google Scholar

Crossref

PubMed

17.

P.

Colli Franzone

,

L. F.

Pavarino

, and

G.

Savaré

, “

Computational electrocardiology: Mathematical and numerical modeling

,” in

Complex Systems in Biomedicine

(

Springer

,

Milan, Milano

,

2006

), pp.

187

–

241

.

Google Scholar

Crossref

18.

P.

Colli Franzone

,

L. F.

Pavarino

, and

S.

Scacchi

, “

Cardiac excitation mechanisms, wavefront dynamics and strength–interval curves predicted by 3d orthotropic bidomain simulations

,”

Math. Biosci.

235

(

1

),

66

–

84

(

2012

).

https://doi.org/10.1016/j.mbs.2011.10.008

Google Scholar

Crossref

PubMed

19.

M. A.

Colman

,

Atrial Fibrillation

(

Springer International Publishing, Cham, Switzerland

,

2014

), pp.

201

–

225

.

Google Scholar

20.

R. L.

DeHaan

and

L. J.

DeFelice

, “

Oscillatory properties and excitability of the heart cell membrane

,”

Theor. Chem.: Period. Chem. Biol.

4

,

181

–

233

(

1978

).

Google Scholar

21.

J.

Engelbrecht

, “

On theory of pulse transmission in a nerve fiber

,”

Proc. R. Soc. London, A

375

,

195

–

209

(

1981

).

https://doi.org/10.1098/rspa.1981.0047

Google Scholar

Crossref

22.

J.

Engelbrecht

,

T.

Peets

,

K.

Tamm

,

M.

Laasmaa

, and

M.

Vendelin

, “

On modelling of physical effects accompanying the propagation of action potentials in nerve fibres

,” preprint arXiv:1601.01867 (

2016

).

Google Scholar

23.

T. E.

Fastl

,

C.

Tobon-Gomez

,

W. A.

Crozier

,

J.

Whitaker

,

R.

Rajani

,

K. P.

McCarthy

,

D.

Sanchez-Quintana

,

S. Y.

Ho

,

M. D.

O’Neill

,

G.

Plank

et al., “

Personalized modeling pipeline for left atrial electromechanics

,” in

2016 Computing in Cardiology Conference (CinC)

(IEEE,

2016

), pp.

225

–

228

.

Google Scholar

Crossref

24.

F. H.

Fenton

,

E. M.

Cherry

,

H. M.

Hastings

, and

S. J.

Evans

, “

Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity

,”

Chaos (Woodbury, N.Y.)

12

(

3

),

852

–

892

(

2002

).

https://doi.org/10.1063/1.1504242

Google Scholar

Crossref

PubMed

25.

F. H.

Fenton

and

A.

Karma

, “

Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation

,”

Chaos

8

(

1

),

20

–

47

(

1998

).

https://doi.org/10.1063/1.166311

Google Scholar

Crossref

PubMed

26.

J.

Fort

,

D.

Campos

,

J. R.

González

, and

J.

Velayos

, “

Bounds for the propagation speed of combustion flames

,”

J. Phys. A: Math. Gen.

37

(

29

),

7185

–

7198

(

2004

).

https://doi.org/10.1088/0305-4470/37/29/002

Google Scholar

Crossref

27.

J.

Fort

and

V.

Méndez

, “

Time-delayed spread of viruses in growing plaques

,”

Phys. Rev. Lett.

89

(

17

),

178101

(

2002a

).

https://doi.org/10.1103/PhysRevLett.89.178101

Google Scholar

Crossref

PubMed

28.

J.

Fort

and

V.

Méndez

, “

Wavefronts in time-delayed reaction-diffusion systems. Theory and comparison to experiment

,”

Rep. Prog. Phys.

65

,

895

(

2002b

).

https://doi.org/10.1088/0034-4885/65/6/201

Google Scholar

Crossref

29.

J.

Gorecki

, “

Molecular dynamics simulations of a chemical wave front

,”

Phys. D: Nonlinear Phenom.

84

(

1–2

),

171

–

179

(

1995

).

https://doi.org/10.1016/0167-2789(95)00033-Z

Google Scholar

Crossref

30.

E.

Grandi

,

S. V.

Pandit

,

N.

Voigt

,

A. J.

Workman

,

D.

Dobrev

,

J.

Jalife

, and

D. M.

Bers

, “

Human atrial action potential and Ca2+ model

,”

Circ. Res.

109

(

9

),

1055

–

1066

(

2011

).

https://doi.org/10.1161/CIRCRESAHA.111.253955

Google Scholar

Crossref

PubMed

31.

B. E.

Griffith

and

C. S.

Peskin

, “

Electrophysiology

,”

Commun. Pure Appl. Math.

66

(

12

),

1837

–

1913

(

2013

).

https://doi.org/10.1002/cpa.21484

Google Scholar

Crossref

PubMed

32.

P. E.

Hand

and

B. E.

Griffith

, “

Adaptive multiscale model for simulating cardiac conduction

,”

Proc. Natl. Acad. Sci.

107

(

33

),

14603

–

14608

(

2010

).

https://doi.org/10.1073/pnas.1008443107

Google Scholar

Crossref

33.

P. E.

Hand

,

B. E.

Griffith

, and

C. S.

Peskin

, “

Deriving macroscopic myocardial conductivities by homogenization of microscopic models

,”

Bull. Math. Biol.

71

(

7

),

1707

–

1726

(

2009

).

https://doi.org/10.1007/s11538-009-9421-y

Google Scholar

Crossref

PubMed

34.

D. M.

Harrild

and

C. S.

Henriquez

, “

A computer model of normal conduction in the human atria

,”

Circ. Res.

87

(

7

),

e25

–

e36

(

2000

).

https://doi.org/10.1161/01.RES.87.7.e25

Google Scholar

Crossref

PubMed

35.

S. Y.

Ho

,

R. H.

Anderson

, and

D.

Sánchez-Quintana

, “

Atrial structure and fibres: Morphologic bases of atrial conduction

,”

Cardiovasc. Res.

54

(

2

),

325

–

336

(

2002

).

https://doi.org/10.1016/S0008-6363(02)00226-2

Google Scholar

Crossref

PubMed

36.

S. Y.

Ho

,

J. A.

Cabrera

, and

D.

Sanchez-Quintana

, “

Left atrial anatomy revisited

,”

Circ.: Arrhythmia Electrophysiol.

5

(

1

),

220

–

228

(

2012

).

https://doi.org/10.1161/CIRCEP.111.962720

Google Scholar

Crossref

37.

A. L.

Hodgkin

and

A. F.

Huxley

, “

A quantitative description of membrane current and its application to conduction and excitation in nerve

,”

J. Physiol.

117

(

4

),

500

(

1952

).

https://doi.org/10.1113/jphysiol.1952.sp004764

Google Scholar

Crossref

PubMed

38.

D. A.

Hooks

,

M. L.

Trew

,

B. J.

Caldwell

,

G. B.

Sands

,

I. J.

LeGrice

, and

B. H.

Smaill

, “

Laminar arrangement of ventricular myocytes influences electrical behavior of the heart

,”

Circ. Res.

101

(

10

),

e103

–

e112

(

2007

).

https://doi.org/10.1161/CIRCRESAHA.107.161075

Google Scholar

Crossref

PubMed

39.

D. E.

Hurtado

,

S.

Castro

, and

A.

Gizzi

, “

Computational modeling of non-linear diffusion in cardiac electrophysiology: A novel porous-medium approach

,”

Comput. Methods Appl. Mech. Eng.

300

,

70

–

83

(

2016

).

https://doi.org/10.1016/j.cma.2015.11.014

Google Scholar

Crossref

40.

D.

Jou

,

J.

Casas-Vázquez

, and

G.

Lebon

, “

Extended irreversible thermodynamics

,” in

Extended Irreversible Thermodynamics

(

Springer

,

Berlin, Heidelberg

,

1996

), pp.

41

–

74

.

Google Scholar

41.

S.

Kaplan

and

D.

Trujillo

, “

Numerical studies of the partial differential equations governing nerve impulse conduction: The effect of Lieberstein’s inductance term

,”

Math. Biosci.

7

,

379

–

404

(

1970

).

https://doi.org/10.1016/0025-5564(70)90135-5

Google Scholar

Crossref

42.

J. P.

Keener

and

J.

Sneyd

,

Mathematical Physiology

(

Springer

,

2009

).

Google Scholar

Crossref

43.

J.

King

,

C.

Huang

, and

J.

Fraser

, “

Determinants of myocardial conduction velocity: Implications for arrhythmogenesis

,”

Front. Physiol.

4

,

154

(

2013

).

https://doi.org/10.3389/fphys.2013.00154

Google Scholar

Crossref

PubMed

44.

B. S.

Kirk

,

J. W.

Peterson

,

R. H.

Stogner

, and

G. F.

Carey

, “

libMesh: A C++ library for parallel adaptive mesh refinement/coarsening simulations

,”

Eng. Comput.

22

(

3–4

),

237

–

254

(

2006

).

https://doi.org/10.1007/s00366-006-0049-3

Google Scholar

Crossref

45.

S. B.

Knisley

, “

Transmembrane voltage changes during unipolar stimulation of rabbit ventricle

,”

Circ. Res.

77

(

6

),

1229

–

1239

(

1995

).

https://doi.org/10.1161/01.RES.77.6.1229

Google Scholar

Crossref

PubMed

46.

C.

Koch

, “

Cable theory in neurons with active, linearized membranes

,”

Biol. Cybern.

50

(

1

),

15

–

33

(

1984

).

https://doi.org/10.1007/BF00317936

Google Scholar

Crossref

PubMed

47.

D.

Krause

,

T.

Dickopf

,

M.

Potse

, and

R.

Krause

, “

Towards a large-scale scalable adaptive heart model using shallow tree meshes

,”

J. Comput. Phys.

298

,

79

–

94

(

2015

).

https://doi.org/10.1016/j.jcp.2015.05.005

Google Scholar

Crossref

48.

S.

Krishnamoorthi

,

L. E.

Perotti

,

N. P.

Borgstrom

,

O. A.

Ajijola

,

A.

Frid

,

A. V.

Ponnaluri

,

J. N.

Weiss

,

Z.

Qu

,

W. S.

Klug

,

D. B.

Ennis

, and

A.

Garfinkel

, “

Simulation methods and validation criteria for modeling cardiac ventricular electrophysiology

,”

PLoS One

9

(

12

),

e114494

(

2014

).

https://doi.org/10.1371/journal.pone.0114494

Google Scholar

Crossref

PubMed

49.

S.

Krishnamoorthi

,

M.

Sarkar

, and

W. S.

Klug

, “

Numerical quadrature and operator splitting in finite element methods for cardiac electrophysiology

,”

Int. J. Numer. Methods Biomed. Eng.

29

(

11

),

1243

–

1266

(

2013

).

https://doi.org/10.1002/cnm.2573

Google Scholar

Crossref

50.

M. W.

Krueger

,

V.

Schmidt

,

C.

Tobón

,

F. M.

Weber

,

C.

Lorenz

,

D. U.

Keller

,

H.

Barschdorf

,

M.

Burdumy

,

P.

Neher

,

G.

Plank

et al., “

Modeling atrial fiber orientation in patient-specific geometries: A semi-automatic rule-based approach

,” in

International Conference on Functional Imaging and Modeling of the Heart

(

Springer

,

2011

), pp.

223

–

232

.

Google Scholar

Crossref

51.

A.

Lemarchand

and

B.

Nawakowski

, “

Perturbation of local equilibrium by a chemical wave front

,”

J. Chem. Phys.

109

(

16

),

7028

–

7037

(

1998

).

https://doi.org/10.1063/1.477685

Google Scholar

Crossref

52.

H. M.

Lieberstein

, “

On the Hodgkin-Huxley partial differential equation

,”

Math. Biosci.

1

(

1

),

45

–

69

(

1967

).

https://doi.org/10.1016/0025-5564(67)90026-0

Google Scholar

Crossref

53.

H. M.

Lieberstein

and

M. A.

Mahrous

, “

A source of large inductance and concentrated moving magnetic fields on axons

,”

Math. Biosci.

7

(

1–2

),

41

–

60

(

1970

).

https://doi.org/10.1016/0025-5564(70)90041-6

Google Scholar

Crossref

54.

MATLAB

,

Version 7.10.0

(

The MathWorks Inc

.,

Natick, Massachusetts

,

2010

).

55.

H. P.

McKean

, “

Nagumo’s equation

,”

Adv. Math.

4

(

3

),

209

–

223

(

1970

).

https://doi.org/10.1016/0001-8708(70)90023-X

Google Scholar

Crossref

56.

V.

Méndez

and

A.

Compte

, “

Wavefronts in bistable hyperbolic reaction-diffusion systems

,”

Physica A

260

,

90

–

98

(

1998

).

https://doi.org/10.1016/S0378-4371(98)00281-7

Google Scholar

Crossref

57.

V.

Méndez

and

J. E.

Llebot

, “

Hyperbolic reaction-diffusion equations for a forest fire model

,”

Phys. Rev. E

56

(

6

),

6557

(

1997

).

https://doi.org/10.1103/PhysRevE.56.6557

Google Scholar

Crossref

58.

M.

Munteanu

,

L. F.

Pavarino

, and

S.

Scacchi

, “

A scalable Newton–Krylov–Schwarz method for the bidomain reaction-diffusion system

,”

SIAM J. Sci. Comput.

31

(

5

),

3861

–

3883

(

2009

).

https://doi.org/10.1137/08074355X

Google Scholar

Crossref

59.

M. P.

Nash

and

A. V.

Panfilov

, “

Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias

,”

Prog. Biophys. Mol. Biol.

85

(

2

),

501

–

522

(

2004

).

https://doi.org/10.1016/j.pbiomolbio.2004.01.016

Google Scholar

Crossref

PubMed

60.

J. C.

Neu

and

W.

Krassowska

, “

Homogenization of syncytial tissues

,”

Crit. Rev. Biomed. Eng.

21

(

2

),

137

–

199

(

1992

).

Google Scholar

61.

S. A.

Niederer

,

E.

Kerfoot

,

A. P.

Benson

,

M. O.

Bernabeu

,

O.

Bernus

,

C.

Bradley

,

E. M.

Cherry

,

R.

Clayton

,

F.

Fenton

,

A.

Garny

,

E.

Heidenreichm

,

S.

Land

,

M.

Maleckar

,

P.

Pathmanathan

,

G.

Plank

,

J. F.

Rodriguez

,

I.

Roy

,

F. B.

Sachse

,

G.

Seemann

,

O.

Skavhaug

, and

N. P.

Smith

, “

Verification of cardiac tissue electrophysiology simulators using an N-version benchmark

,”

Philos. Trans.. Ser. A, Math., Phys. Eng. Sci.

369

,

4331

–

4351

(

2011

).

https://doi.org/10.1098/rsta.2011.0139

Google Scholar

Crossref

62.

J. W.

Papez

, “

Heart musculature of the atria

,”

Dev. Dyn.

27

(

3

),

255

–

285

(

1920

).

https://doi.org/10.1002/aja.1000270302

Google Scholar

Crossref

63.

F.

Pashakhanloo

,

D. A.

Herzka

,

H.

Ashikaga

,

S.

Mori

,

N.

Gai

,

D. A.

Bluemke

,

N. A.

Trayanova

, and

E. R.

McVeigh

, “

Myofiber architecture of the human atria as revealed by submillimeter diffusion tensor imaging

,”

Circ.: Arrhythmia Electrophysiol.

9

(

4

),

e004133

(

2016

).

https://doi.org/10.1161/CIRCEP.116.004133

Google Scholar

Crossref

64.

A. S.

Patelli

,

L.

Dede

,

T.

Lassila

,

A.

Bartezzaghi

, and

A.

Quarteroni

, “

Isogeometric approximation of cardiac electrophysiology models on surfaces: An accuracy study with application to the human left atrium

,”

Comput. Methods Appl. Mech. Eng.

317

,

248

–

273

(

2017

).

https://doi.org/10.1016/j.cma.2016.12.022

Google Scholar

Crossref

65.

P.

Pathmanathan

,

G. R.

Mirams

,

J.

Southern

, and

J. P.

Whiteley

, “

The significant effect of the choice of ionic current integration method in cardiac electro-physiological simulations

,”

Int. J. Numer. Methods Biomed. Eng.

27

(

11

),

1751

–

1770

(

2011

).

https://doi.org/10.1002/cnm.1438

Google Scholar

Crossref

66.

S.

Pezzuto

,

J.

Hake

, and

J.

Sundnes

, “

Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology

,”

Int. J. Numer. Methods Biomed. Eng.

32

(

10

),

e02762

(

2016

).

https://doi.org/10.1002/cnm.2762

Google Scholar

Crossref

67.

M.

Potse

,

B.

Dube

,

J.

Richer

,

A.

Vinet

, and

R. M.

Gulrajani

, “

A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart

,”

IEEE Trans. Biomed. Eng.

53

(

12

),

2425

–

2435

(

2006

).

https://doi.org/10.1109/TBME.2006.880875

Google Scholar

Crossref

PubMed

68.

A.

Quarteroni

,

T.

Lassila

,

S.

Rossi

, and

R.

Ruiz-Baier

, “

Integrated heart-coupling multiscale and multiphysics models for the simulation of the cardiac function

,”

Comput. Methods Appl. Mech. Eng.

314

,

345

–

407

(

2017

).

https://doi.org/10.1016/j.cma.2016.05.031

Google Scholar

Crossref

69.

S.

Rossi

,

T.

Lassila

,

R.

Ruiz-Baier

,

A.

Sequeira

, and

A.

Quarteroni

, “

Thermodynamically consistent orthotropic activation model capturing ventricular systolic wall thickening in cardiac electromechanics

,”

Eur. J. Mech. A/Solids

48

,

129

–

142

(

2014

).

https://doi.org/10.1016/j.euromechsol.2013.10.009

Google Scholar

Crossref

70.

B. J.

Roth

, “

Electrical conductivity values used with the bidomain model of cardiac tissue

,”

IEEE Trans. Biomed. Eng.

44

(

4

),

326

–

328

(

1997

).

https://doi.org/10.1109/10.563303

Google Scholar

Crossref

PubMed

71.

S.

Rush

and

H.

Larsen

, “

A practical algorithm for solving dynamic membrane equations

,”

IEEE Trans. Biomed. Eng.

BME-25

(

4

),

389

–

392

(

1978

).

https://doi.org/10.1109/TBME.1978.326270

Google Scholar

Crossref

72.

F. B.

Sachse

,

Computational Cardiology

, Lecture Notes in Computer Science (

Springer

,

Berlin, Heidelberg

,

2004

), Vol.

2966

.

Google Scholar

Crossref

73.

D.

Sánchez-Quintana

,

G.

Pizarro

,

J. R.

López-Mínguez

,

S. Y.

Ho

, and

J. A.

Cabrera

, “

Standardized review of atrial anatomy for cardiac electrophysiologists

,”

J. Cardiovasc. Transl. Res.

6

(

2

),

124

–

144

(

2013

).

https://doi.org/10.1007/s12265-013-9447-2

Google Scholar

Crossref

PubMed

74.

A. C.

Scott

, “

Effect of the series inductance of a nerve axon upon its conduction velocity

,”

Math. Biosci.

11

(

3–4

),

277

–

290

(

1971

).

https://doi.org/10.1016/0025-5564(71)90089-7

Google Scholar

Crossref

75.

A. C.

Scott

, “

Transmission line equivalent for an unmyelinated nerve axon

,”

Math. Biosci.

13

(

1–2

),

47

–

54

(

1972

).

https://doi.org/10.1016/0025-5564(72)90022-3

Google Scholar

Crossref

76.

W. P.

Segars

,

G.

Sturgeon

,

S.

Mendonca

,

J.

Grimes

, and

B. M. W.

Tsui

, “

4d xcat phantom for multimodality imaging research

,”

Med. Phys.

37

(

9

),

4902

–

4915

(

2010

).

https://doi.org/10.1118/1.3480985

Google Scholar

Crossref

PubMed

77.

N. G.

Sepulveda

,

B. J.

Roth

, and

J. P.

Wikswo

, “

Current injection into a two-dimensional anisotropic bidomain

,”

Biophys. J.

55

(

5

),

987

–

999

(

1989

).

https://doi.org/10.1016/S0006-3495(89)82897-8

Google Scholar

Crossref

PubMed

78.

M. S.

Spach

, “

The discontinuous nature of electrical propagation in cardiac muscle—Consideration of a quantitative model incorporating the membrane ionic properties and structural complexities the ALZA distinguished lecture

,”

Ann. Biomed. Eng.

11

(

3–4

),

208

–

261

(

1983

).

https://doi.org/10.1007/BF02363287

Google Scholar

Crossref

79.

R. J.

Spiteri

and

R. C.

Dean

, “

Stiffness analysis of cardiac electrophysiological models

,”

Ann. Biomed. Eng.

38

(

12

),

3592

–

3604

(

2010

).

https://doi.org/10.1007/s10439-010-0100-9

Google Scholar

Crossref

PubMed

80.

D.

Stan

,

F.

del Teso

, and

J. L.

Vázquez

, “

Finite and infinite speed of propagation for porous medium equations with fractional pressure

,”

C. R. Math.

352

(

2

),

123

–

128

(

2014

).

https://doi.org/10.1016/j.crma.2013.12.003

Google Scholar

Crossref

81.

K. H. W. J.

ten Tusscher

and

A. V.

Panfilov

, “

Alternans and spiral breakup in a human ventricular tissue model

,”

Am. J. Physiol.: Heart Circ. Physiol.

291

(

3

),

H1088

–

H1100

(

2006

).

Google Scholar

PubMed

82.

C. E.

Thomas

, “

The muscular architecture of the atria of hog and dog hearts

,”

Dev. Dyn.

104

(

2

),

207

–

236

(

1959

).

Google Scholar

83.

L.

Tung

,

A Bi-Domain Model for Describing Ischemic Myocardial D-C Potentials

(

Diss. Massachusetts Institute of Technology

,

1978

).

Google Scholar

84.

Z.

Vera

,

H. P.

Pride

,

D. P.

Zipes

, and

R. C.

Barr

, “

Reperfusion arrhythmias: Role of early afterdepolarizations studied by monophasic action potential recordings in the intact canine heart during autonomically denervated and stimulated states

,”

J. Cardiovasc. Electrophysiol.

6

(

7

),

532

–

543

(

1995

).

https://doi.org/10.1111/j.1540-8167.1995.tb00426.x

Google Scholar

Crossref

PubMed

85.

J. P.

Wikswo

,

S. F.

Lin

, and

R. A.

Abbas

, “

Virtual electrodes in cardiac tissue: A common mechanism for anodal and cathodal stimulation

,”

Biophys. J.

69

(

6

),

2195

–

2210

(

1995

).

https://doi.org/10.1016/S0006-3495(95)80115-3

Google Scholar

Crossref

PubMed

86.

J. P.

Wikswo

and

B. J.

Roth

, “

Virtual electrode theory of pacing

,” in

Cardiac Bioelectric Therapy

(

Springer

,

2009

), pp.

283

–

330

.

Google Scholar

Crossref

87.

E. P.

Zemskov

,

M. A.

Tsyganov

, and

W.

Horsthemke

, “

Wavy fronts in a hyperbolic FitzHugh-Nagumo system and the effects of cross diffusion

,”

Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.

91

(

6

),

062917

(

2015

).

https://doi.org/10.1103/PhysRevE.91.062917

Google Scholar

Crossref

2017

Author(s)

Incorporating inductances in tissue-scale models of cardiac electrophysiology

I. PHENOMENOLOGICAL DERIVATION OF THE HYPERBOLIC BIDOMAIN EQUATIONS

II. REDUCTION TO THE HYPERBOLIC MONODOMAIN

III. IONIC MODELS

IV. NUMERICAL RESULTS

A. Comparison with an exact solution

B. Effect of the relaxation time on the conduction velocity for different ionic models

C. Conduction velocity anisotropy

D. Discretization error: A simple two-dimensional benchmark

E. Effect of the relaxation time on spiral break up

F. Atrial fibrillation

G. Virtual electrodes

V. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: CABLE EQUATION WITH INDUCTANCES

APPENDIX B: EXACT SOLUTION FOR THE PIECEWISE-LINEAR BISTABLE MODEL

APPENDIX C: FINITE ELEMENT DISCRETIZATION

References

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Incorporating inductances in tissue-scale models of cardiac electrophysiology Free

I. PHENOMENOLOGICAL DERIVATION OF THE HYPERBOLIC BIDOMAIN EQUATIONS

II. REDUCTION TO THE HYPERBOLIC MONODOMAIN

III. IONIC MODELS

IV. NUMERICAL RESULTS

A. Comparison with an exact solution

B. Effect of the relaxation time on the conduction velocity for different ionic models

C. Conduction velocity anisotropy

D. Discretization error: A simple two-dimensional benchmark

E. Effect of the relaxation time on spiral break up

F. Atrial fibrillation

G. Virtual electrodes

V. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: CABLE EQUATION WITH INDUCTANCES

APPENDIX B: EXACT SOLUTION FOR THE PIECEWISE-LINEAR BISTABLE MODEL

APPENDIX C: FINITE ELEMENT DISCRETIZATION

References

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Incorporating inductances in tissue-scale models of cardiac electrophysiology