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,37 the cable equation and its higher-dimensional generalizations83 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,53 and hyperbolic models including inductances have been proposed by Engelbrecht,21 Lieberstein,52 and Engelbrecht et al.22 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 Defelice15 and Koch.46 The common conclusion that inductances are negligible, which is based on the linear analyses for neurons by Scott74,75 and the numerical results of Kaplan and Trujillo,41 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 Lieberstein52 and Lieberstein and Mahrous53 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,41 Scott,74 and Engelbrecht,21 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.75 Several years later, further studies by DeHaan and DeFelice,20 Clapham and Defelice,15 and Koch46 showed that the cell membranes of excitable tissue exhibit self-inductance. In particular, Clapham and Defelice15 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.87 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 Cattaneo14 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,57 chemical systems,29,51 thermal combustion,26 and the spread of viral infections.27 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.48 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,66 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.33 As shown by Quarteroni et al.,68 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.32 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.

In the 1970s, Tung83 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.,13 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.31 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, {f0,s0,n0}, we assume the conductivity tensors can be written as σj=σjff0f0+σjss0s0+σjnn0n0, for j=e,i. Typical values of the conductivity coefficients range from 102 mS/cm in the cross-fiber direction to 1 mS/cm in the fiber direction (for example, Colli Franzone et al.17 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

τeJet+Je=σeVe,
(1)
τiJit+Ji=σiVi,
(2)

where Vi (Ve), Ji (Je), 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.,40 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 equation28,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.

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).

FIG. 1.

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).

Close modal

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

τet·Je+·Je=·σeVe,
(3)
τit·Ji+·Ji=·σiVi.
(4)

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

·(Ji+Je)=0.
(5)

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

·Ji=It=·Je,
(6)

where It=χ(CmVt+Iion) is the usual transmembrane current density, with χ being the membrane area per unit volume of tissue, Cm the membrane capacitance, and Iion the transmembrane ionic current. Using Eq. (6) in (3) and (4), we obtain the system of equations

τeItt+It=·σeVe,
(7)
τiIttIt=·σiVi.
(8)

Defining the transmembrane potential V=ViVe to eliminate Vi from the equations yields

τeItt+It=·σeVe,
(9)
τiItt+It=·σiV+·σiVe.
(10)

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

τeCm2Vt2+CmVt+·DeVe=IionτeIiont,
(11)
τiCm2Vt2+CmVt·DiV·DiVe=IionτiIiont,
(12)

where we set Di=σi/χ and De=σe/χ. Alternatively, the model equations can be written as

τiCm2Vt2+CmVt·DiV·DiVe=IionτiIiont,
(13)
(τeτi)Cm2Vt2+·DiV+·(De+Di)Ve=(τiτe)Iiont.
(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

τCm2Vt2+CmVt·DiV·DiVe=IionτIiont,
(15)
·DiV+·(De+Di)Ve=0.
(16)

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

CmVt·DiV·DiVe=Iion,
(17)
·DiV+·(De+Di)Ve=0.
(18)

The hyperbolic bidomain model can be simplified by assuming that the extracellular and intracellular compartments have the same anisotropy ratios, so that D=De=λDi. If we make this assumption, we obtain from (14)

·DVe=1λ+1·DV+λλ+1(τeτi)Cm2Vt2λλ+1(τiτe)Iiont.
(19)

Substituting in (13), we find

[τi+λλ+1(τeτi)]Cm2Vt2+CmVtλλ+1·DV=Iion[τi+λλ+1(τeτi)]Iiont.
(20)

Defining τ=τi+λ(τeτi)/(λ+1), we obtain the hyperbolic monodomain model

τCm2Vt2+CmVt·DV=IionτIiont,
(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=Vt, we transform the hyperbolic monodomain equations into the first-order system

Vt=Q,
(22)
τCmQt+CmQ·DV=IionτIiont.
(23)

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

V·N=0
(24)

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

Following Stan et al.,80 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(·,t) is also compactly supported for any t > 0. The compact support is taken with respect to the resting potential V0. 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 MR,t={x:||x||2R,V(x,t)>V0} has a positive measure. For any relaxation time τ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.

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 80mV. The bidomain and monodomain models must be completed by specifying the form of Iion, 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 Iion=Iion(V,w). Typically, the state variable dynamics are determined by a spatially decoupled system of nonlinear ordinary differential equations

wt=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,55 the two-variable model of Aliev and Panfilov (AP),3 the three-variable model of Fenton and Karma,25 and the biophysically detailed models of Tusscher and Panfilov81 for the ventricles (20 variables) and of Grandi et al.30 for the atria (57 variables).

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

rt=(ε+μ1rμ2+V)(rkV(Vb1)),
(26)

and the ionic current is Iion=kV(Vα)(V1)+rV. The parameters used in this case are shown in Table I. In the three-variable model of Fenton and Karma,25 w={v,w}, with

vt=1τv(V)(1p)(1v)1τv+pv,
(27)
wt=1τw(1p)(1w)1τw+pw,
(28)

where τv(V)=(1q)τv1+qτv2,p=H(VVc),q=H(VVv), and H(·) is the Heaviside function. The total ionic current is the sum of three currents, Iion=Ifi(V,v)Iso(V)Isi(V,w), with

Ifi(V,v)=1τdvp(VVc)(1V),
(29)
Iso(V)=1τ0V(1p)+1τrp,
(30)
Isi(V,w)=12τsiw(1+tanh(k(VVcsi))).
(31)

The sets of parameters used in this paper for this model are shown in Table II. We refer to the original papers30,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).

Cmσfσs=σnχkbμ1μ2ε
0.125 0.1 0.12 0.3 0.01 
Cmσfσs=σnχkbμ1μ2ε
0.125 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.24 

Cmσfσs=σnχτv+τv1τv2τw+τwτdτ0τrτsikVcsiVcVv
Parameter set 3 0.1 0.0125 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 0.1 0.0125 3.33 15.6 350 80 0.407 9.0 34.0 26.5 15 0.45 0.15 0.04 
Parameter set 5 0.1 0.0125 3.33 12.0 1000 100 0.362 5.0 33.33 29.0 15 0.7 0.13 0.04 
Parameter set 6 0.1 0.0125 3.33 9.0 250 60 0.395 9.0 33.33 29.0 15 0.5 0.13 0.04 
Cmσfσs=σnχτv+τv1τv2τw+τwτdτ0τrτsikVcsiVcVv
Parameter set 3 0.1 0.0125 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 0.1 0.0125 3.33 15.6 350 80 0.407 9.0 34.0 26.5 15 0.45 0.15 0.04 
Parameter set 5 0.1 0.0125 3.33 12.0 1000 100 0.362 5.0 33.33 29.0 15 0.7 0.13 0.04 
Parameter set 6 0.1 0.0125 3.33 9.0 250 60 0.395 9.0 33.33 29.0 15 0.5 0.13 0.04 

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 library44 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 Matlab54 to visualize the one-dimensional results and Paraview2 for the two- and three-dimensional simulations.

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

Iion(V)=kVk[V2H(VV1)+V0(1H(VV1))],
(32)

where V0 is the resting potential, V1 is the threshold potential, V2 is the depolarization potential, and H(·) is the Heaviside function. This model (32) reduces to the piecewise-linear model of McKean55 after a simple dimensional analysis, after which the nondimensional ionic currents take the form

Îion(V̂)=V̂H(V̂α).
(33)

The nondimensional potential V̂ is related to the dimensional potential by V=(V2V0)V̂+V0. In this model, α=(V1V0)/(V2V0) 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=(12α)μ+(αα2)(μ1)2<1μ=cs.
(34)

The speed cs=1/μ=Cm/τ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/Cm 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=σkχCm2(12α)μ+(αα2)(μ1)2<σχτCm.
(35)

As a verification test, we consider the spatial interval Ω=[0,50]. An initial stimulus of amplitude 1 is applied at x[24.5,25.5] for the interval t[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 V̂ and Q̂ at t = 10 in Fig. 2 (center and right). We register the activation time at a particular spatial location whenever V̂ crosses a threshold of 0.9. The conduction velocities are measured by picking the distance between two points x1 and x2Ω 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 x1=30 and x2=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.

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

FIG. 2.

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

Close modal

We also examine the effect of the relaxation time on the conduction velocities for three values of the excitability parameter, α{0.1,0.2,0.3}. With this simple ionic model, the larger the relaxation time, the slower the wave speed. The limit speed cs 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

Îion(V̂)=V̂(V̂1)(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.

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.

FIG. 3.

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.

Close modal

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,3 the three-variable model of Fenton and Karma,25 the ventricular M-cell model of Tusscher and Panfilov81 (20 variables), and the atrial model of Grandi et al.30 (57 variables). For the two-variable Aliev–Panfilov model, we use the parameters of Nash and Panfilov,59 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.24 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 Cm=1μF/cm2 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[2.45,2.55] cm, during the time interval t[0.03,1.03] ms. The conduction velocities are measured, as in the previous test case, by dividing the distance between two points x1 and x2 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.

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 τ.

FIG. 4.

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 τ.

Close modal

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.

2DMax CG iterationsReaction step costDiffusion step costOverall
100 × 100τ = 0 msτ = 0.4 msτ=0.4ms/τ=0 msτ=0.4ms/τ=0 msτ=0.4ms/τ=0 ms
AP −54% 23% −3% 
FK 22 15 −62% 22% −2% 
TP06 11 −7% 15% −3% 
2DMax CG iterationsReaction step costDiffusion step costOverall
100 × 100τ = 0 msτ = 0.4 msτ=0.4ms/τ=0 msτ=0.4ms/τ=0 msτ=0.4ms/τ=0 ms
AP −54% 23% −3% 
FK 22 15 −62% 22% −2% 
TP06 11 −7% 15% −3% 
3DMax CG iterationsReaction step costDiffusion step costOverall
10 × 10×10τ = 0 msτ = 0.4 msτ=0.4ms/τ=0τ=0.4ms/τ=0 msτ=0.4ms/τ=0 ms
AP 43 16 −73% 48% 30% 
FK 43 16 −49% 47% 25% 
TP06 43 16 −6% 46% 15% 
3DMax CG iterationsReaction step costDiffusion step costOverall
10 × 10×10τ = 0 msτ = 0.4 msτ=0.4ms/τ=0τ=0.4ms/τ=0 msτ=0.4ms/τ=0 ms
AP 43 16 −73% 48% 30% 
FK 43 16 −49% 47% 25% 
TP06 43 16 −6% 46% 15% 

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 v1, v2, v3, and v4 corresponding to the conductivities σ1, σ2, σ3, and σ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 II, σ1 and σ4 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 σ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 τ=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 (σ1:σ4), 4:1 (σ1:σ3), and 2:1 (σ1:σ2), we find that the ratios between the conduction velocities in the transverse and in the longitudinal directions are approximately 1/8 (v4/v1), 1/2 (v3/v1), and 1/2 (v2/v1). 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.

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.

FIG. 5.

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.

Close modal

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]×[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Ω:||x||11} 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 h234μm and h117μm show the greatest differences for transverse propagation.

FIG. 6.

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. 6.

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.

Close modal
FIG. 7.

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.

FIG. 7.

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.

Close modal

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.,68 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.32 

Our next tests explore the effect of relaxation time on spiral wave break up. We consider a square slab of tissue, Ω=[0,12]×[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<6cmy<7cm} 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.

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]×[0,12] cm, the mesh size is h100μm, and the time step size Δt=0.125 ms. The fibers are aligned with the y-axis.

FIG. 8.

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]×[0,12] cm, the mesh size is h100μm, and the time step size Δt=0.125 ms. The fibers are aligned with the y-axis.

Close modal

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=tV. This variable can be used to define the fronts (red) and the tails (light blue) of the waves.

FIG. 9.

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]×[0,12] cm, the mesh size is h200μm in the fiber direction, and h50μm in the transverse direction, and the time step is Δt=0.125 ms. The fibers are aligned with the y-axis.

FIG. 9.

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]×[0,12] cm, the mesh size is h200μm in the fiber direction, and h50μm in the transverse direction, and the time step is Δt=0.125 ms. The fibers are aligned with the y-axis.

Close modal

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.67 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.76 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 × 106 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 Thomas82) 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.69 In the other regions, similar to work by Patelli et al.64 and Krueger et al.,50 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.

(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.

FIG. 10.

(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.

Close modal

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.19 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,34 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.

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.

FIG. 11.

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.

Close modal

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 Ie,stim in the extracellular compartment, the region of hyperpolarization near a cathode (Ie,stim<0) is called a virtual anode, while the region of depolarization near an anode (Ie,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 Roth86 and Colli Franzone et al.18 

Following Sepulveda et al.,77 we consider the domain Ω=[2,2]cm×[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×[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 σef=1.5448, σes=σen=1.0438, and σif=2.3172, σis=σin=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.

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.

FIG. 12.

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.

Close modal

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.39 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.39 use the simplified ionic model proposed by Aliev and Panfilov.3 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.,43 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,55 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 myocytes81 and atrial myocytes30 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.25 

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.,24 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.,35 Krueger et al.,50 and Pashakhanloo et al.63 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 velocity48,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.68 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 literature49,65 and have been thoroughly analyzed by Pezzuto et al.,66 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.61 Other solution strategies for the monodomain and bidomain models using adaptive mesh refinement17,47 and high-order elements4 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.

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.

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

Ve(x)Ve(x+Δx)=LetIeΔxReIeΔx,
(A1)
Vi(x)Vi(x+Δx)=LitIiΔxRiIiΔx,
(A2)

where Ie and Ii are the extracellular and intracellular axial currents, respectively. Dividing by Δx and taking the limit as Δx0, we find

xVe=LetIeReIe,
(A3)
xVi=LitIiRiIi.
(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

Ie(x)Ie(x+Δx)+ItΔx=0,
(A5)
Ii(x)Ie(x+Δx)ItΔx=0,
(A6)

which, in the limit Δx0, becomes

It=xIe=xIi.
(A7)

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

xxVe=LextIeRexIe,
(A8)
xtIe=tIt,
(A9)

and for the intracellular compartment,

xxVi=LixtIiRixIi,
(A10)
xtIe=tIt.
(A11)

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

xxVe=LetItReIt,xxVi=LitIt+RiIt.
(A12)

Defining V=ViVe, and taking the difference between the intracellular and the extracellular equations, we can write a single equation for V

xxV=(Li+Le)tIt+(Ri+Re)It.
(A13)

Recalling that

It=χ(CmtV+Iion),
(A14)

we finally find the hyperbolic form of the cable equation

τCm2Vt2+CmVtx(DVx)=IionτIiont,

where we define the conductivity as D=1/χ(Ri+Re) and the relaxation time by

τ=Li+LeRi+Re.
(A15)

Consider the simplified piecewise-linear bistable model

Iion(V)=kVk[V2H(VV1)+V0(1H(VV1))].
(B1)

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

Îion(V̂)=V̂H(V̂α).
(B2)

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

{μ2V̂t̂2+(1+μ)V̂t̂2V̂x̂2+V̂=0,V̂<α,μ2V̂t̂2+(1+μ)V̂t̂2V̂x̂2+V̂=1,V̂>α,
(B3)

where we have introduced the nondimensional variables

t̂=1Tt̂,x̂=1Lx,V̂=VV0V2V0.

In particular, we define

T=Cmk,L=σkχ,μ=τkCm,

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=xct, such that U(z)=U(xct)=V̂(x,t), the system (B3) is transformed into

{(c2μ1)Uzzc(1+μ)Uz+U=0,z>0,(c2μ1)Uzzc(1+μ)Uz+U=1,z<0.
(B4)

Consider first the case z > 0. Using γ=c2μ1 and β=c(1+μ), this equation reads

γV+βV+V=0

for which the roots of the characteristic polynomial are

λ±=β±β24γ2γ,
(B5)

so that the solution is U(z)=A+eλ+z+Aeλz. It is easy to verify that the case z < 0 has the solution U(z)=B+eλ+z+Beλz+1. The global solution, therefore, is

U(z)={A+eλ+z+Aeλz,z>0,B+eλ+z+Beλz+1,z<0.
(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 β24γ>0, we find the condition γ<β2/4, which is always satisfied for the model (B4) because β24γ=c2(1μ)2+4>0. If c > 0 then β<0. Then λ>0 and λ+<0 if

c<1μ.
(B7)

Imposing U()=0, we find A=0, and imposing U()=1, we find B+=0, so that

U(z)={A+eλ+z,z>0,Beλz+1,z0.
(B8)

To ensure the continuity of the solution at z = 0, we fix U(0)=α, so that A+=α=1+B and

U(z)={αeλ+z,z>0,(α1)eλz+1,z0.
(B9)

The first derivative of U is

U(z)={αλ+eλ+z,z>0,(α1)λeλz,z0.
(B10)

Imposing the continuity of U at z = 0, we find αλ+=(α1)λ and therefore

γβ2=(α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=(12α)μ+(αα2)(μ1)21μ.
(B12)

Notice that the characteristic propagation speed is bounded from above by 1/μ.

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.,58 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 [tn,tn+1], with t0=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 wn+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,79 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,81 we use the Rush–Larsen method.71 The atrial ionic model of Grandi et al.30 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=tn+1tn. 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)Ω=Ωvw the L2(Ω) inner product and using the boundary conditions (24), the Galerkin approximation of Eqs. (22) and (23) is to find Vh,QhSh such that

(Vht,ψh)Ω=(Qh,ψh)Ω,
(C1)
(τCmQht+CmQh,ϕh)Ω+(DV,ϕh)Ω=(Iionh,ϕh)Ω+(τIionht,ϕh)Ω
(C2)

for all ψh,ϕhSh, where

Sh={vhC0(Ω¯):vh|KP1(K),KTh}.
(C3)

In practice, we look for a continuous solution that is linear in every simplex element K in the triangulation Th of Ω. Introducing the basis of nodal shape functions {NA}A=1M, with M=dim(Sh), the solution fields are discretized as

Vh=A=1MNAVA,Qh=A=1MNAQA.

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

V̇=Q,
(C4)
CmM(τQ̇+Q)+KV=F+τL.
(C5)

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

F=MI,L=MJ

and

V̇=Q,
(C6)
CmML(τQ̇+Q)+KV=M(I+τJ),
(C7)

where

[MAB]=(NA,NB)Ω,[KAB]=(NA,DNB)Ω,

I and J are the nodal evaluation of Iion and Iion/t, and ML 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

Vn+1=Vn+ΔtQn+1,
(C8)
Cm(τ+Δt)MLQn+1+ΔtKVn+1=ΔtM(I*+τJ*)+τCmMLQn,
(C9)

where I* and J* indicate the values of the ionic currents, and of their derivatives, evaluated using the updated values for wn+1 obtained by solving the ionic model. Introducing (C8) in (C9), we arrive at a single system for Qn+1,

[Cm(τ+Δt)ML+Δt2K]Qn+1=τCmMLQnΔtKVn+ΔtM(I*+τJ*).
(C10)

The time derivative of the ionic currents, J*, is approximated nodally using the chain rule as

Iiont(tn)IionV(Vn,wn+1)Qn+Iionw(Vn,wn+1)·g(Vn,wn+1),
(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,10 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 hm the smallest element size, we chose our time step size such that the CFL condition Δthm/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.

1.
H. A.
Abdusalam
, “
Analytic and approximate solutions for Nagumo telegraph reaction diffusion equation
,”
Appl. Math. Comput.
157
(
2
),
515
522
(
2004
).
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
.
3.
R. R.
Aliev
and
A. V.
Panfilov
, “
A simple two-variable model of cardiac excitation
,”
Chaos, Solitons Fractals
7
(
3
),
293
301
(
1996
).
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
).
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
).
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
).
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
.
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
.
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
).
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
).
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
).
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
).
14.
C.
Cattaneo
, “
Sulla conduzione del calore
,”
Atti Semin. Mat. Fis. Univ. Modena
3
,
83
101
(
1948
).
15.
D. E.
Clapham
and
L. J.
Defelice
, “
Membrane biolow small signal impedance of heart cell membranes
,”
J. Membr. Biol.
67
,
63
71
(
1982
).
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
).
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
.
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
).
19.
M. A.
Colman
,
Atrial Fibrillation
(
Springer International Publishing, Cham, Switzerland
,
2014
), pp.
201
225
.
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
).
21.
J.
Engelbrecht
, “
On theory of pulse transmission in a nerve fiber
,”
Proc. R. Soc. London, A
375
,
195
209
(
1981
).
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
).
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
.
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
).
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
).
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
).
27.
J.
Fort
and
V.
Méndez
, “
Time-delayed spread of viruses in growing plaques
,”
Phys. Rev. Lett.
89
(
17
),
178101
(
2002a
).
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
).
29.
J.
Gorecki
, “
Molecular dynamics simulations of a chemical wave front
,”
Phys. D: Nonlinear Phenom.
84
(
1–2
),
171
179
(
1995
).
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
).
31.
B. E.
Griffith
and
C. S.
Peskin
, “
Electrophysiology
,”
Commun. Pure Appl. Math.
66
(
12
),
1837
1913
(
2013
).
32.
P. E.
Hand
and
B. E.
Griffith
, “
Adaptive multiscale model for simulating cardiac conduction
,”
Proc. Natl. Acad. Sci.
107
(
33
),
14603
14608
(
2010
).
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
).
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
).
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
).
36.
S. Y.
Ho
,
J. A.
Cabrera
, and
D.
Sanchez-Quintana
, “
Left atrial anatomy revisited
,”
Circ.: Arrhythmia Electrophysiol.
5
(
1
),
220
228
(
2012
).
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
).
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
).
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
).
40.
D.
Jou
,
J.
Casas-Vázquez
, and
G.
Lebon
, “
Extended irreversible thermodynamics
,” in
Extended Irreversible Thermodynamics
(
Springer
,
Berlin, Heidelberg
,
1996
), pp.
41
74
.
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
).
42.
J. P.
Keener
and
J.
Sneyd
,
Mathematical Physiology
(
Springer
,
2009
).
43.
J.
King
,
C.
Huang
, and
J.
Fraser
, “
Determinants of myocardial conduction velocity: Implications for arrhythmogenesis
,”
Front. Physiol.
4
,
154
(
2013
).
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
).
45.
S. B.
Knisley
, “
Transmembrane voltage changes during unipolar stimulation of rabbit ventricle
,”
Circ. Res.
77
(
6
),
1229
1239
(
1995
).
46.
C.
Koch
, “
Cable theory in neurons with active, linearized membranes
,”
Biol. Cybern.
50
(
1
),
15
33
(
1984
).
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
).
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
).
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
).
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
.
51.
A.
Lemarchand
and
B.
Nawakowski
, “
Perturbation of local equilibrium by a chemical wave front
,”
J. Chem. Phys.
109
(
16
),
7028
7037
(
1998
).
52.
H. M.
Lieberstein
, “
On the Hodgkin-Huxley partial differential equation
,”
Math. Biosci.
1
(
1
),
45
69
(
1967
).
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
).
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
).
56.
V.
Méndez
and
A.
Compte
, “
Wavefronts in bistable hyperbolic reaction-diffusion systems
,”
Physica A
260
,
90
98
(
1998
).
57.
V.
Méndez
and
J. E.
Llebot
, “
Hyperbolic reaction-diffusion equations for a forest fire model
,”
Phys. Rev. E
56
(
6
),
6557
(
1997
).
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
).
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
).
60.
J. C.
Neu
and
W.
Krassowska
, “
Homogenization of syncytial tissues
,”
Crit. Rev. Biomed. Eng.
21
(
2
),
137
199
(
1992
).
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
).
62.
J. W.
Papez
, “
Heart musculature of the atria
,”
Dev. Dyn.
27
(
3
),
255
285
(
1920
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
70.
B. J.
Roth
, “
Electrical conductivity values used with the bidomain model of cardiac tissue
,”
IEEE Trans. Biomed. Eng.
44
(
4
),
326
328
(
1997
).
71.
S.
Rush
and
H.
Larsen
, “
A practical algorithm for solving dynamic membrane equations
,”
IEEE Trans. Biomed. Eng.
BME-25
(
4
),
389
392
(
1978
).
72.
F. B.
Sachse
,
Computational Cardiology
, Lecture Notes in Computer Science (
Springer
,
Berlin, Heidelberg
,
2004
), Vol.
2966
.
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
).
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
).
75.
A. C.
Scott
, “
Transmission line equivalent for an unmyelinated nerve axon
,”
Math. Biosci.
13
(
1–2
),
47
54
(
1972
).
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
).
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
).
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
).
79.
R. J.
Spiteri
and
R. C.
Dean
, “
Stiffness analysis of cardiac electrophysiological models
,”
Ann. Biomed. Eng.
38
(
12
),
3592
3604
(
2010
).
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
).
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
).
82.
C. E.
Thomas
, “
The muscular architecture of the atria of hog and dog hearts
,”
Dev. Dyn.
104
(
2
),
207
236
(
1959
).
83.
L.
Tung
,
A Bi-Domain Model for Describing Ischemic Myocardial D-C Potentials
(
Diss. Massachusetts Institute of Technology
,
1978
).
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
).
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
).
86.
J. P.
Wikswo
and
B. J.
Roth
, “
Virtual electrode theory of pacing
,” in
Cardiac Bioelectric Therapy
(
Springer
,
2009
), pp.
283
330
.
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
).