Incorporating Inductances in Tissue-Scale Models of Cardiac Electrophysiology

TL;DR

This paper introduces a hyperbolic bidomain and monodomain model for cardiac electrophysiology that incorporates inductances, addressing the unrealistic infinite propagation speed in standard models and exploring nonlinear effects on conduction velocity.

Contribution

It develops a hyperbolic formulation of cardiac electrophysiology models that includes inductances, providing a more realistic representation of wave propagation and nonlinear dynamics.

Findings

Inductances can increase conduction velocity in nonlinear models.

Hyperbolic models require finer spatial discretization for accurate velocity prediction.

The models influence spiral wave dynamics and atrial fibrillation simulations.

Abstract

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 [...].

Tables

Table 1. Table 1: Parameters for the Aliev-Panfilov ionic model ( 26 ).

$C_{\text{m}}$	${\sigma }^{\text{f}}$	${\sigma }^{\text{s}}={\sigma }^{\text{n}}$	$\chi$	$k$	$b$	${\mu }_1$	${\mu }_2$	$\epsilon$
1	1	0.125	1	8	0.1	0.12	0.3	0.01

Table 2. Table 2: Sets of parameters for the Fenton-Karma ionic model ( 27 )–( 31 ) as stated by Fenton et al. ( 2002 ) .

	$C_{\text{m}}$	${\sigma }^{\text{f}}$	${\sigma }^{\text{s}}={\sigma }^{\text{n}}$	$\chi$	${\tau }_v^{+}$	${\tau }_{v1}^{-}$	${\tau }_{v2}^{-}$	${\tau }_w^{+}$	${\tau }_w^{-}$	${\tau }_{\text{d}}$	${\tau }_0$	${\tau }_{\text{r}}$	${\tau }_{\text{si}}$	$k$	$V_{\text{c}}^{\text{si}}$	$V_{\text{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

Table 3. Table 3: Comparison of the compuational cost for the parabolic ( τ = 0 𝜏 0 \tau=0 ms) and hyperbolic ( τ = 0.4 𝜏 0.4 \tau=0.4 ms) monodomain models for two- and three-dimensional problems in a serial computation. Negative values indicate relative slowdown of the hyperbolic model compared to the parabolic model, whereas 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$\times$100	$\tau =0$ ms	$\tau =0$.4 ms	$\tau =0.4\text{ ms}/\tau =0$ ms	$\tau =0.4\text{ ms}/\tau =0$ ms	$\tau =0.4\text{ ms}/\tau =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$\times$10$\times$10	$\tau =0$ ms	$\tau =0$.4 ms	$\tau =0.4\text{ ms}/\tau =0$	$\tau =0.4\text{ ms}/\tau =0$ ms	$\tau =0.4\text{ ms}/\tau =0$ ms
AP	43	16	-73%	48%	30%
FK	43	16	-49%	47%	25%
TP06	43	16	-6%	46%	15%

Equations

$${\tau }_{\text{e}}\frac{\partial {\mathbf{J}}_{\text{e}}}{\partial t}+{\mathbf{J}}_{\text{e}}$$

$${\tau }_{\text{i}}\frac{\partial {\mathbf{J}}_{\text{i}}}{\partial t}+{\mathbf{J}}_{\text{i}}$$

$${\tau }_{\text{e}}\frac{\partial }{\partial t}\nabla \cdot {\mathbf{J}}_{\text{e}}+\nabla \cdot {\mathbf{J}}_{\text{e}}$$

$${\tau }_{\text{i}}\frac{\partial }{\partial t}\nabla \cdot {\mathbf{J}}_{\text{i}}+\nabla \cdot {\mathbf{J}}_{\text{i}}$$

$$\nabla \cdot \left({\mathbf{J}}_{\text{i}}+{\mathbf{J}}_{\text{e}}\right)=0.$$

$$-\nabla \cdot {\mathbf{J}}_{\text{i}}=I_{\text{t}}=\nabla \cdot {\mathbf{J}}_{\text{e}},$$

$${\tau }_{\text{e}}\frac{\partial I_{\text{t}}}{\partial t}+I_{\text{t}}$$

$$-{\tau }_{\text{i}}\frac{\partial I_{\text{t}}}{\partial t}-I_{\text{t}}$$

$${\tau }_{\text{e}}\frac{\partial I_{\text{t}}}{\partial t}+I_{\text{t}}$$

$${\tau }_{\text{i}}\frac{\partial I_{\text{t}}}{\partial t}+I_{\text{t}}$$

$${\tau }_{\text{e}}{C}_{\text{m}}\frac{{\partial }^{2}V}{\partial t^2}+C_{\text{m}}\frac{\partial V}{\partial t}+\nabla \cdot {\mathbf{D}}_{\text{e}}\nabla V_{\text{e}}=$$

$$-I_{\text{ion}}-{\tau }_{\text{e}}\frac{\partial I_{\text{ion}}}{\partial t},$$

$${\tau }_{\text{i}}{C}_{\text{m}}\frac{{\partial }^{2}V}{\partial t^2}+C_{\text{m}}\frac{\partial V}{\partial t}$$

$$-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V_{\text{e}}=$$

$$-I_{\text{ion}}-{\tau }_{\text{i}}\frac{\partial I_{\text{ion}}}{\partial t},$$

$${\tau }_{\text{i}}{C}_{\text{m}}\frac{{\partial }^{2}V}{\partial t^2}+C_{\text{m}}\frac{\partial V}{\partial t}$$

$$-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V_{\text{e}}=$$

$$-I_{\text{ion}}-{\tau }_{\text{i}}\frac{\partial I_{\text{ion}}}{\partial t},$$

$$\left({\tau }_{\text{e}}-{\tau }_{\text{i}}\right)C_{\text{m}}\frac{{\partial }^{2}V}{\partial t^2}$$

$$+\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V+\nabla \cdot \left({\mathbf{D}}_{\text{e}}+{\mathbf{D}}_{\text{i}}\right)\nabla V_{\text{e}}=$$

$$\left({\tau }_{\text{i}}-{\tau }_{\text{e}}\right)\frac{\partial I_{\text{ion}}}{\partial t}.$$

$$\tau C_{\text{m}}\frac{{\partial }^{2}V}{\partial t^2}+C_{\text{m}}\frac{\partial V}{\partial t}$$

$$-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V_{\text{e}}$$

$$\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V+\nabla \cdot \left({\mathbf{D}}_{\text{e}}+{\mathbf{D}}_{\text{i}}\right)\nabla V_{\text{e}}$$

$$C_{\text{m}}\frac{\partial V}{\partial t}-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V-\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V_{\text{e}}$$

$$\nabla \cdot {\mathbf{D}}_{\text{i}}\nabla V+\nabla \cdot \left({\mathbf{D}}_{\text{e}}+{\mathbf{D}}_{\text{i}}\right)\nabla V_{\text{e}}$$

$$-\nabla \cdot \mathbf{D}\nabla V_{\text{e}}=\frac{1}{\lambda +1}\nabla \cdot \mathbf{D}\nabla V$$

$$+\frac{\lambda }{\lambda +1}\left({\tau }_{\text{e}}-{\tau }_{\text{i}}\right)C_{\text{m}}\frac{{\partial }^{2}V}{\partial t^2}-\frac{\lambda }{\lambda +1}\left({\tau }_{\text{i}}-{\tau }_{\text{e}}\right)\frac{\partial I_{\text{ion}}}{\partial t},$$

$$\left[{\tau }_{\text{i}}+\frac{\lambda }{\lambda +1}\left({\tau }_{\text{e}}-{\tau }_{\text{i}}\right)\right]C_{\text{m}}\frac{{\partial }^{2}V}{\partial t^2}$$

$$+C_{\text{m}}\frac{\partial V}{\partial t}-\frac{\lambda }{\lambda +1}\nabla \cdot \mathbf{D}\nabla V=-I_{\text{ion}}$$

$$-\left[{\tau }_{\text{i}}+\frac{\lambda }{\lambda +1}\left({\tau }_{\text{e}}-{\tau }_{\text{i}}\right)\right]\frac{\partial I_{\text{ion}}}{\partial t}.$$

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

$$\frac{\partial V}{\partial t}$$

$$\tau C_{\text{m}}\frac{\partial Q}{\partial t}+C_{\text{m}}Q-\nabla \cdot \mathbf{D}\nabla V$$

$$\nabla V\cdot \mathbf{N}=0$$

$$\frac{\partial \mathbf{w}}{\partial t}=\mathbf{g}\left(V,\mathbf{w}\right),$$

$$\frac{\partial r}{\partial t}=\left(\epsilon +\frac{{\mu }_{1}r}{{\mu }_2+V}\right)\left(-r-kV\left(V-b-1\right)\right),$$

$$\frac{\partial v}{\partial t}$$

$$\frac{\partial w}{\partial t}$$

$$I_{\text{fi}}\left(V,v\right)$$

$$I_{\text{so}}\left(V\right)$$

Full text

Cardiac Electrophysiology, Spiral waves, Finite Element Method

Since its introduction by Hodgkin & Huxley (1952), the cable equation and its higher-dimensional generalizations Tung (1978) 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 & Mahrous (1970), and hyperbolic models including inductances have been proposed by Lieberstein (1967); Engelbrecht (1981), and Engelbrecht et al. (2016). These models are also supported by the fact that neurons, skeletal muscle cells, and cardiomyocytes show the typical resonance effects due to inductances, as demonstrated by the studies of Clapham & Defelice (1982) and Koch (1984). The common conclusion that inductances are negligible, which is based on the linear analyses for neurons by Scott Scott (1971, 1972) and the numerical results of Kaplan & Trujillo (1970), 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 cell 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 Spach (1983); Vera et al. (1995); Bueno-Orovio et al. (2008). 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 works of Lieberstein (1967) and Lieberstein & Mahrous (1970) 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 & Trujillo (1970); Scott (1971), and Engelbrecht (1981), have investigated this hypothesis for nerves. Their conclusions were that inductances in neurons are of the order of few $\mu$H and therefore are negligible, following a linear analysis of a one-dimensional nerve discussed by Scott (1972). Several years later, further studies by DeHaan & DeFelice (1978); Clapham & Defelice (1982), and Koch (1984) showed that the cell membranes of excitable tissue exhibit self-inductance. In particular, Clapham & Defelice (1982) 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. (2015) 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 (1948) 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 Méndez & Llebot (1997), chemical systems Gorecki (1995); Lemarchand & Nawakowski (1998), thermal combustion Fort et al. (2004), and the spread of viral infections Fort & Méndez (2002a). 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 Jou et al. (1996); Méndez & Compte (1998). 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. (2014) propose that the “wave speed should not be sensitive to choices of 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 Pezzuto et al. (2016), the common belief within the field seems to be that spatial discretizations on the order of 200 $\mu$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 Hand et al. (2009). As shown by Quarteroni et al. (2017), the mesh sizes needed to resolve transverse propagation are actually closer to 25 $\mu$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 Hand & Griffith (2010). 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 1970’s, Tung (1978) 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. The bidomain equations can be derived from a model that accounts for the tissue microarchitecture Keener & Sneyd (2009); Neu & Krassowska (1992); Hand et al. (2009), but 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. (2014), 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 & Peskin (2013). 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 $\boldsymbol{\sigma}_{\text{e}}$ and $\boldsymbol{\sigma}_{\text{i}}$. Introducing a local orthonormal basis, $\left\{\boldsymbol{f}_{\text{0}},\boldsymbol{s}_{\text{0}},\boldsymbol{n}_{\text{0}}\right\}$, we assume the conductivity tensors can be written as $\boldsymbol{\sigma}_{j}=\sigma_{j}^{\text{f}}\boldsymbol{f}_{\text{0}}\otimes\boldsymbol{f}_{\text{0}}+\sigma_{j}^{\text{s}}\boldsymbol{s}_{\text{0}}\otimes\boldsymbol{s}_{\text{0}}+\sigma_{j}^{\text{n}}\boldsymbol{n}_{\text{0}}\otimes\boldsymbol{n}_{\text{0}},$ for $j=\mathrm{e},\,\mathrm{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. (2006) use $\sigma^{\text{f}}=1.2$ mS/cm, $\sigma^{\text{s}}=0.346$ mS/cm, and $\sigma^{\text{n}}=0.0435$ mS/cm). We model the current density fluxes using a Cattaneo-type equation, so that

[TABLE]

where $V_{\text{i}}$ ($V_{\text{e}}$), $\boldsymbol{J}_{\text{i}}$ ($\boldsymbol{J}_{\text{e}}$), and $\tau_{\text{i}}$ ($\tau_{\text{e}}$) are the intracellular (extracellular) potential, the intracellular (extracellular) flux, and the intracellular (extracellular) relaxation time, respectively. As shown by Jou et al. (1996), 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 Jou et al. (1996); Fort & Méndez (2002b) shows that they are consistent with the theory of extended irreversible thermodynamics. Alternatively, equations (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 equation (50) the relationship between inductance and the relaxation time.

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

[TABLE]

As in similar derivations of the bidomain model, we impose a quasistatic form of charge conservation Sachse (2004), yielding

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where we set $\boldsymbol{D}_{\text{i}}=\boldsymbol{\sigma}_{\text{i}}/\chi$ and $\boldsymbol{D}_{\text{e}}=\boldsymbol{\sigma}_{\text{e}}/\chi$. Alternatively, the model equations can be written as

[TABLE]

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

[TABLE]

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

[TABLE]

II Reduction to the hyperbolic monodomain

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

[TABLE]

Substituting in (13), we find

[TABLE]

Defining $\tau=\tau_{\text{i}}+\lambda\left(\tau_{\text{e}}-\tau_{\text{i}}\right)/\left(\lambda+1\right)$, we obtain the hyperbolic monodomain model,

[TABLE]

where here we have absorbed the term $\lambda/\left(\lambda+1\right)$ into $\boldsymbol{D}$. The relaxation time $\tau$ of the monodomain model is always positive, and it is zero only if both $\tau_{\text{i}}$ and $\tau_{\text{e}}$ are zero. In fact, if $\tau_{\text{e}}=0$, then $\tau=\tau_{\text{i}}/\left(\lambda+1\right)$, and if $\tau_{\text{i}}=0$, then $\tau=\lambda\tau_{\text{e}}/\left(\lambda+1\right)$. However, if $\tau_{\text{i}}=\tau_{\text{e}}$, then $\tau=\tau_{\text{i}}=\tau_{\text{e}}$.

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

[TABLE]

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

[TABLE]

on the boundary of the domain, where the vector $\boldsymbol{N}$ is the normal to the boundary.

Following Stan et al. (2014), 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_{0}$. 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}=\{\boldsymbol{x}:||\boldsymbol{x}||_{2}\geq R,V(\boldsymbol{x},t)>V_{0}\}$ has positive measure. For any relaxation time $\tau\neq 0$, (21) is hyperbolic, and it thereby has the property of finite propagation speed. Setting $\tau=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\leavevmode\nobreak\ \text{mV}$. The bidomain and monodomain models must be completed by specifying the form of $I_{\text{ion}}$, which accounts for these transmembrane currents. The states of the transmembrane ion channels are described by a collection of variables, $\boldsymbol{w}$, associated with the ionic model, so that $I_{\text{ion}}=I_{\text{ion}}(V,\boldsymbol{w})$. Typically, the state variable dynamics are determined by a spatially decoupled system of nonlinear ordinary differential equations,

[TABLE]

where the form of $\boldsymbol{g}$ depends on the details of the particular ionic model. We consider the simplified piecewise-linear model of McKean (1970), the two-variable model of Aliev & Panfilov (1996), the three-variable model of Fenton & Karma (1998), and the biophysically detailed models of ten Tusscher & Panfilov (2006) for the ventricles (20 variables) and of Grandi et al. (2011) for the atria (57 variables).

In the two-variable model of Aliev & Panfilov (1996), $\boldsymbol{w}=\left\{r\right\}$, with

[TABLE]

and the ionic current is $I_{\text{ion}}=kV\left(V-\alpha\right)\left(V-1\right)+rV$. The parameters used in this case are shown in Table 1. In the three-variable model of Fenton & Karma (1998), $\boldsymbol{w}=\left\{v,w\right\}$, with

[TABLE]

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

[TABLE]

The sets of parameters used in this paper for this model are shown in Table 2. We refer to the original papers ten Tusscher & Panfilov (2006); Grandi et al. (2011) for the statements of the equations and parameters of the biophysically detailed ionic models.

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(Kirk et al., 2006) and relies on linear solvers provided by PETSc (Balay et al., 1997, 2015a, 2015b). All the code used for the following tests is 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(MATLAB, 2010) to visualize the one-dimensional results and Paraview(Ahrens et al., 2005) for the two- and three-dimensional simulations.

IV.1 Comparison with an exact solution

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

[TABLE]

where $V_{0}$ is the resting potential, $V_{1}$ is the threshold potential, $V_{2}$ is the depolarization potential, and $H\left(\cdot\right)$ is the Heaviside function. This model (32) reduces to the piecewise-linear model of McKean (1970) after a simple dimensional analysis, after which the nondimensional ionic currents take the form

[TABLE]

The nondimensional potential $\hat{V}$ is related to the dimensional potential by $V=\left(V_{2}-V_{0}\right)\hat{V}+V_{0}$. In this model, $\alpha=\left(V_{1}-V_{0}\right)/\left(V_{2}-V_{0}\right)$ 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

[TABLE]

The speed $c_{\text{s}}=\sqrt{1/\mu}=\sqrt{C_{\text{m}}/\tau k}$ represents the maximum speed at which a perturbation can travel in the system. When the relaxation time $\tau$ goes to zero, the local perturbations can travel at infinite speed. The parameter $\mu=\tau k/C_{\text{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

[TABLE]

As a verification test, we consider the spatial interval $\Omega=\left[0,50\right]$. An initial stimulus of amplitude $1$ is applied at $x\in\left[24.5,25.5\right]$ for the interval $t\in\left[0.03,1.03\right]$. The system of equations is solved using a mesh size $h=0.03125$ and a time step size $\Delta 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_{1}$ and $x_{2}$ $\in\Omega$ and dividing it by the time interval between the activation times at these two locations. Because (34) is obtained by assuming that $\Omega$ is the entire real line, we need to measure 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).

We also examine the effect of the relaxation time on the conduction velocities for three values of the excitability parameter, $\alpha\in\left\{0.1,\,0.2,\,0.3\right\}$. With this simple ionic model, the larger the relaxation time, the slower the wave speed. The limit speed $c_{\text{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 $\Omega=[-25,25]$ and set the initial conditions according to (59) and (60), assuming the front at time $t=0$ is at $x=0$ and $\alpha=0.1$. The system is discretized in time using implicit-explicit (IMEX) Runge-Kutta (RK) time integratorsAscher et al. (1997); Boscarino et al. (2016). 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 (59) and (60) 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 $\Delta t$ is taken to be linearly proportional to the mesh size $h$, with a value of $\Delta t=0.05$ in the coarsest cases. The time step size was chosen to guarantee a large enough number of time iterations and while satisfying the CFL condition (see Appendix C). Fig. 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 derivative 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-$\delta$ functions. The errors reported in Fig. 3 correspond to the more 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

[TABLE]

The exact solution for the cubic reaction in the parabolic monodomain has been determined previously Keener & Sneyd (2009); Pezzuto et al. (2016). 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.

IV.2 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 & Panfilov (1996), the three-variable model of Fenton & Karma (1998), the ventricular M-cell model of ten Tusscher & Panfilov (2006) (20 variables), and the atrial model of Grandi et al. (2011) (57 variables). For the two-variable Aliev-Panfilov model, we use the parameters of Nash & Panfilov (2004), as reported in Table 1. We consider two values for the excitability parameter, $\alpha=0.1$ and $\alpha=0.2$. For the three-variable Fenton-Karma model, we consider parameters sets 3, 4, 5, and 6 given by Fenton et al. (2002) and reported in Table 2. 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_{\text{m}}=1$ $\mu\text{F}/\text{cm}^{2}$ for the membrane capacitance in both models.

For all cases, the domain is the interval $\Omega=\left[0,5\right]$ cm. The simulations use a mesh size $h=31.25$ $\mu$m and a time step size $\Delta t=0.003653$ ms. An initial stimulus is applied at the center of the domain, $x\in\left[2.45,2.55\right]$ cm, during the time interval $t\in\left[0.03,1.03\right]$ ms. The conduction velocities are measured, as in the previous test case, by dividing the distance between two points $x_{1}$ and $x_{2}$ 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 $\left[0,1\right]$ 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.

We also compare the computational time required by the parabolic and hyperbolic monodomain models in two and three spatial dimensions. Table 3 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 3 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 system 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$.

IV.3 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: $\sigma_{1}=0.1$ mS/mm, $\sigma_{2}=0.05$ mS/mm, $\sigma_{3}=0.025$ mS/mm, and $\sigma_{4}=0.0125$ mS/mm. Evaluating the wave speeds $v_{1}$, $v_{2}$, $v_{3}$, and $v_{4}$ corresponding to the conductivities $\sigma_{1}$, $\sigma_{2}$, $\sigma_{3}$, and $\sigma_{4}$, respectively, we show how the anisotropy ratio changes with respect to the relaxation time for several conductivity ratios. In fact, as shown in Table 2, $\sigma_{1}$ and $\sigma_{4}$ are the typical longitudinal and transverse conductivities for this model. Considering the interval $\Omega=\left[0,5\right]$ cm, we use a mesh size $h=25$ $\mu$m and a time step size $\Delta t=0.0025$ ms. Fig. 5(left) shows how the velocities at different $\sigma$ 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 $\tau=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. Fig. 5(right) shows the relaxation time needed to maintain the same conduction velocity as in the parabolic monodomain model. In particular, for $\sigma_{4}$, 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 $\tau=0.38$ ms. The differences in the velocities computed with $\tau=0.4$ ms and $\tau=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 $\tau=0.4$ ms. Fig. 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 ($\sigma_{1}$:$\sigma_{4}$), 4:1 ($\sigma_{1}$:$\sigma_{3}$), and 2:1 ($\sigma_{1}$:$\sigma_{2}$), we find that the ratios between the conduction velocities in the transverse and in the longitudinal directions are approximately 1/$\sqrt{8}$ ($v_{4}$/$v_{1}$), 1/$2$ ($v_{3}$/$v_{1}$), and 1/$\sqrt{2}$ ($v_{2}$/$v_{1}$). 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.

IV.4 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, $\Omega=\left[0,12\right]\times\left[0,12\right]$ 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, $\Omega_{\text{stim}}=\left\{\boldsymbol{x}\in\Omega:\left\|\boldsymbol{x}\right\|_{1}\leq 1\right\}$ 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. Figs. 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, $\Omega_{\text{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 $\mu$m (163 $\mu$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 $\mu$m (82.5 $\mu$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\text{234 }\mu\text{m}$ and $h\approx\text{117 }\mu\text{m}$ show the greatest differences for transverse propagation.

Although similar results have already been reported in prior work Krishnamoorthi et al. (2013); Niederer et al. (2011); Pezzuto et al. (2016), it is nonetheless generally believed that a mesh size of the order of 200 $\mu$m is usually sufficient for cardiac computational electrophysiologyKrishnamoorthi et al. (2014); Fastl et al. (2016). This belief is supported by numerical simulations using a mesh size of 200 $\mu$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. (2017), the required grid resolution to capture the transverse conduction velocity with an error smaller than 5% is about 25 $\mu$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 & Griffith (2010).

IV.5 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, $\Omega=\left[0,12\right]\times\left[0,12\right]$ 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 $\left\{x<6\text{ cm}\,\wedge\,y<7\text{ cm}\right\}$ 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 2, 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.

Fig. 8 shows the formation of the spiral wave in the case of isotropy. With $\tau=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. 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 $\tau=0.4$ ms. As explained in Sec. IV.3, although the conduction velocities in the transverse direction are different in the hyperbolic and parabolic models for $\tau=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. Fig. 9 shows the dynamics of $Q=\partial_{\text{t}}V$. This variable can be used to define the fronts (red) and the tails (light blue) of the waves.

IV.6 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 equations (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. (2006) 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. (2010) 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 million 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 Papez (1920); Pashakhanloo et al. (2016); Sánchez-Quintana et al. (2013), but the structure of the muscle in the atria is so complex (as can be seen from the diagrams by Thomas (1959)) 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. (2014) In the other regions, similar to work by Patelli et al. (2017) and Krueger et al. (2011), 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. Ho et al. (2002, 2012). 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 direction.

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 (2014). 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. Fig. 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 Harrild & Henriquez (2000), 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 $\tau=0$ ms only after the third stimulus is applied. Fig. 11 also shows the fronts (red) and tails (light blue) of the waves, as indicated by $Q$.

IV.7 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 studiesKnisley (1995); Wikswo et al. (1995), 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_{\text{e,stim}}$ in the extracellular compartment, the region of hyperpolarization near a cathode ($I_{\text{e,stim}}<0$) is called a virtual anode, while the region of depolarization near an anode ($I_{\text{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 & Roth (2009) and Colli Franzone et al. (2012).

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

Fig. 12 shows the pattern formed by the transmembrane potential 2 ms after the extracellular stimulus is applied. Red regions are areas of depolarization, and blue regions are 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 $\tau_{\text{i}}$ and $\tau_{\text{e}}$ to take the values 0, 0.2, and 0.4 ms. Note that if $\tau_{\text{i}}=\tau_{\text{e}}$, the hyperbolic bidomain equations (13) and (14) reduce to the hyperbolic-elliptic system of equations (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.

V Conclusions

Local perturbation in the standard bidomain and monodomain models propagate 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 equations (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. (2016) 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. (2016) use the simplified ionic model proposed by Aliev & Panfilov (1996). 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. (2013), 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 (1970), 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 Abdusalam (2004); Méndez & Compte (1998).) 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 Aliev & Panfilov (1996); Fenton & Karma (1998), 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 ten Tusscher & Panfilov (2006) and atrial myocytes Grandi et al. (2011) 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 Fenton & Karma (1998).

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 do 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 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 conduction block, leading to a break up of the wave. This behavior, already observed for the standard monodomain model and explained by Fenton et al. (2002), 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. (2002); Krueger et al. (2011), and Pashakhanloo et al. (2016) 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 $\mu$m is sufficient to correctly capture the conduction velocity Krishnamoorthi et al. (2013, 2014) 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 $\mu$m (163 $\mu$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 $\mu$m (82.5 $\mu$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 Roth (1997); Clayton et al. (2011); Hooks et al. (2007); Benson et al. (2010), the mesh size needs to be much smaller than 200 $\mu$m in order to yield resolved dynamics. Our results are in accordance with the convergence study by Rossi et al. (2014) for propagation in the transverse direction, in which the authors showed that the mesh size in the transverse direction needs to be about 25 $\mu$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 Pathmanathan et al. (2011); Krishnamoorthi et al. (2013) and have been thoroughly analyzed by Pezzuto et al. (2016), 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. (2011) Other solution strategies for the monodomain and bidomain models using adaptive mesh refinement Colli Franzone et al. (2006); Krause et al. (2015) and high-order elements Arthurs et al. (2012) 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.

Acknowledgements

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

Appendix A Cable equation with inductances

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

[TABLE]

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

[TABLE]

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 that

[TABLE]

which, in the limit $\Delta x\to 0$, becomes

[TABLE]

Differentiating equations (38) and (39) with respect to $x$, and equation (42) with respect to $t$, we find that for the extracellular compartment,

[TABLE]

and for the intracellular compartment,

[TABLE]

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

[TABLE]

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

[TABLE]

Recalling that

[TABLE]

we finally find the hyperbolic form of the cable equation,

[TABLE]

where we define the conductivity as $D=1/\chi\left(R_{\text{i}}+R_{\text{e}}\right)$ and the relaxation time by

[TABLE]

Appendix B Exact solution for the piecewise-linear bistable model

Consider the simplified piecewise-linear bistable model,

[TABLE]

The nondimensional form of equation (51) is the simplified model proposed by McKean (1970),

[TABLE]

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

[TABLE]

where we have introduced the nondimensional variables

[TABLE]

In particular, we define

[TABLE]

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

[TABLE]

Consider first the case $z>0$. Using $\gamma=c^{2}\mu-1$ and $\beta=-c\left(1+\mu\right)$, this equation reads

[TABLE]

for which the roots of the characteristic polynomial are

[TABLE]

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

[TABLE]

For (56) to represent a propagating front, it is necessary that the solution is real and bounded for any value of $z$. This implies that either $\lambda_{-}<0<\lambda_{+}$ or $\lambda_{+}<0<\lambda_{-}$. Requiring $\beta^{2}-4\gamma>0$, we find the condition $\gamma<\beta^{2}/4$, which is always satisfied for the model (54) because $\beta^{2}-4\gamma=c^{2}\left(1-\mu\right)^{2}+4>0.$ If $c>0$ then $\beta<0$. Then $\lambda_{-}>0$ and $\lambda_{+}<0$ if

[TABLE]

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

[TABLE]

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

[TABLE]

The first derivative of $U$ is

[TABLE]

Imposing the continuity of $U^{\prime}$ at $z=0$, we find $\alpha\lambda_{+}=\left(\alpha-1\right)\lambda_{-}$ and therefore

[TABLE]

Using the definitions of $\gamma$ and $\beta$ in (61), we finally obtain an expression for the speed $c$ in terms of $\mu$ and $\alpha$,

[TABLE]

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

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. (2009), 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 $\boldsymbol{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 & Dean (2010), the most efficient numerical algorithm to solve the stiff ODEs 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 ten Tusscher & Panfilov (2006), we use the Rush-Larsen method Rush & Larsen (1978). The atrial ionic model of Grandi et al. (2011) 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 $\Delta 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)–(23) is comparable. Denoting by $\left(v,w\right)_{\Omega}=\int_{\Omega}vw$ the $L^{2}\left(\Omega\right)$ inner product and using the boundary conditions (24), the Galerkin approximation of equations (22)–(23) is to find $V^{h},\,Q^{h}$ $\in\mathcal{S}^{h}$ such that

[TABLE]

for all $\psi^{h},$ $\phi^{h}$ $\in\mathcal{S}^{h}$, where

[TABLE]

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

[TABLE]

Equations (63) and (64) yield the matrix system

[TABLE]

We use the half-lumping scheme proposed by Pathmanathan et al. (2011), so that

[TABLE]

and

[TABLE]

where

[TABLE]

$\boldsymbol{\mathsf{I}}$ and $\boldsymbol{\mathsf{J}}$ are the nodal evaluation of $I_{\text{ion}}$ and $\partial I_{\text{ion}}/\partial t$, and $\boldsymbol{\mathsf{M}}_{\text{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

[TABLE]

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

[TABLE]

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

[TABLE]

where the function $\boldsymbol{g}(V,\boldsymbol{w})$ is the right hand side of the ionic model system defined in equation (25). After solving (72), we update the voltage equation (70).

The scheme used to obtain (70) and (71) is the simplest method – ARS(1,1,1) – of a series of implicit-explicit (IMEX) Runge-Kutta (RK) algorithms that are popular for hyperbolic systems Ascher et al. (1997); Boscarino & Russo (2009). For a second-order time integrator, we use the H-CN(2,2,2) schemeBoscarino et al. (2016), 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 an 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 $\Delta 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.