Incorporating Inductances in Tissue-Scale Models of Cardiac Electrophysiology
IIncorporating Inductances in Tissue-Scale Models of CardiacElectrophysiology
Simone Rossi a) and Boyce E. Griffith Department of Mathematics, University of North Carolina, Chapel Hill, NC USA Departments of Mathematics and Biomedical Engineering and McAllister Heart Institute,University of North Carolina, Chapel Hill, NC USA
In standard models of cardiac electrophysiology, including the bidomain and monodomain models, local per-turbations can propagate at infinite speed. We address this unrealistic property by developing a hyperbolicbidomain 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 con-ductivity tensors have the same anisotropy ratio. In one spatial dimension, the hyperbolic monodomain modelis equivalent to a cable model that includes axial inductances, and the relaxation times of the Cattaneo fluxesare 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 fullynonlinear 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 sim-ilar behavior is also demonstrated with biophysically detailed ionic models. Using the Fenton-Karma modelalong with a low-order finite element spatial discretization, we numerically analyze differences between thestandard monodomain model and the hyperbolic monodomain model. In a simple benchmark test, we showthat the propagation of the action potential is strongly influenced by the alignment of the fibers with respectto 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 singlecardiac cell. We also compare the two formulations in the case of spiral break up and atrial fibrillation in ananatomically detailed model of the left atrium, and we examine the effect of intracellular and extracellularinductances on the virtual electrode phenomenon.Keywords: Cardiac Electrophysiology, Spiral waves, Finite Element Method
Since its introduction by Hodgkin and Huxley ,the cable equation and its higher-dimensionalgeneralizations have been common models ofelectrical impulse propagation in excitable me-dia, including neurons and muscle. The effectsof inductances in these systems are consideredto be relatively small, and so they are neglectedin classical versions of these models. By omit-ting such terms, the standard equations of car-diac electrophysiology become parabolic, but, asin all parabolic equations, local perturbationscan propagate at infinite speed. This unreal-istic property has been addressed in models ofneurons by Lieberstein and Mahrous , and hy-perbolic models including inductances have beenproposed by Engelbrecht , Lieberstein , andEngelbrecht et al. . These models are also sup-ported by the fact that neurons, skeletal mus-cle cells, and cardiomyocytes show the typicalresonance effects due to inductances, as demon-strated by the studies of Clapham and Defelice and Koch . The common conclusion that induc-tances are negligible, which is based on the lin-ear analyses for neurons by Scott and the nu-merical results of Kaplan and Trujillo , may not a) Email: [email protected] be valid in the complex arrangement of cardiactissue, where the inhomogeneities together withhighly nonlinear reactions can lead to reentrantwaves and chaotic behavior. We thus derive a hy-perbolic model for cardiac electrophysiology, andwe compare the solutions with parabolic modelsin several cases, including simple excitation pat-terns as well as for spiral waves, atrial fibrillation,and the virtual electrode phenomenon.
The heart is a complex organ with a highly heteroge-neous structure. Muscle cell are embedded in an extracel-lular compartment that includes many components, in-cluding capillaries, connective tissue, and collagen. Thestructural arrangement of the tissue is known to influ-ence electrical impulse propagation . The stan-dard models of electrophysiology neglect all these com-plexities, and the resulting equations have the unrealisticfeature that local perturbations can propagate infinitelyfast. These models are fundamentally based on the ca-ble equation. The works of Lieberstein and Liebersteinand Mahrous were the first to suggest that the cablemodel for neurons should contain inductances because ofthe three-dimensional nature of the axon. Several au-thors, including Kaplan and Trujillo , Scott , and En-gelbrecht , have investigated this hypothesis for nerves.Their conclusions were that inductances in neurons areof the order of few µ H and therefore are negligible, fol-lowing a linear analysis of a one-dimensional nerve dis- a r X i v : . [ m a t h . NA ] A ug cussed by Scott . Several years later, further studies byClapham and Defelice , DeHaan and DeFelice , andKoch showed that the cell membranes of excitable tis-sue exhibit self-inductance. In particular, Clapham andDefelice showed that the impedance of the embryonicheart cell membrane resonates at a frequency around 1Hz, thereby enhancing homogeneity of the voltage. Toour knowledge, there has been no experimental studythat addresses the role of inductances for reentrant waves.In this chaotic scenario, inductances may have an impor-tant role in the system dynamics.Excitable tissues have a characteristic speed of trans-mission that limits the velocity at which signals can prop-agate. Any signal, including action potential wavefronts,cannot propagate faster than this characteristic speed.On the other hand, it is important to notice that propa-gation of the action potential is related to the nonlineardynamics of the system, and not to the wave-like behaviorthat may be induced by inductances. One way to trans-form the cable equation into a hyperbolic equation is tosimply add an “inertial” term proportional to the persis-tence time of the diffusive process, as done by Zemskovet al. Unfortunately, this type of model cannot be de-rived from reasonable physical assumptions. However, aswe show here, it is possible to perform a phenomenologi-cal derivation of a hyperbolic reaction-diffusion model byusing the Cattaneo model for the fluxes. This formula-tion was originally introduced by Cattaneo to elim-inate the anomalies found using Fourier’s law in theheat equation, and this model has subsequently beenused in a wide range of applications, including forest firemodels , chemical systems , thermal combustion ,and the spread of viral infections . Cattaneo-type mod-els for the fluxes have been derived in several ways, rang-ing from phenomenological and thermodynamical deriva-tions to isotropic and anisotropic random walks withreactions . The use of Cattaneo-type fluxes in a mon-odomain model of cardiac electrophysiology leads to twoadditional terms in the equations proportional to thecharacteristic relaxation time of the medium: the sec-ond derivative in time of the potential, which is associ-ated with “inertia”, and the time derivative of the ioniccurrents. Even if the relaxation time is small, the rapidvariation of the ionic currents can impart a contributionthat may not be negligible. This is particularly relevantnear the front of the wave, where fast currents give riseto the upstroke of the action potential.Verification and validation remain major challengesin computational electrophysiology. Krishnamoorthiet al. propose that the “wave speed should not be sen-sitive to choices of numerical solution protocol, such asmesh density, numerical integration scheme, etc.” Al-though this is a fundamentally desirable criterion, it isalso nearly impossible to achieve in practice. What doesseem reasonable is to require that that the error in thewave speed be “small”, but how small the error mustbe taken clearly depends on the case under considera-tion. For example, for a single heart beat, an error of 5% might represent a reasonable approximation. On theother hand, models of cardiac fibrillation may require amuch smaller error. In fact, in this case, a 5% error inthe conduction velocity can determine whether reentrantwaves form, or when and how they break up. Althoughthis problem has been discussed in detail in prior work ,the common belief within the field seems to be that spa-tial discretizations on the order of 200 µ m are sufficientto capture the conduction velocities of the propagatingfronts. This estimate seems to hold for isotropic propaga-tion at normal coupling strengths, but in many importantcases, the most relevant propagation of the electrical sig-nal is transverse to the alignment of the cardiac cells,where the conductivities are typically 8 times smallerthan in the longitudinal direction . As shown by Quar-teroni et al. , the mesh sizes needed to resolve trans-verse propagation are actually closer to 25 µ m. Such asmall mesh spacing requires the use of large-scale simula-tions and highly efficient codes. Further, the need to usesuch high spatial resolutions challenges the fundamen-tal idea of using a continuum model for the descriptionof cardiac electrophysiology, and multiscale models havebeen proposed . This paper shows that the transverseconduction velocities are sensitive to the grid size andto the mesh orientation for both regular and hyperbolicversions of the model, and that the hyperbolic model hassimilar mesh size requirements as the standard model. I. PHENOMENOLOGICAL DERIVATION OF THEHYPERBOLIC BIDOMAIN EQUATIONS
In the 1970’s, Tung formulated a bidomain modelof the propagation of the action potential in cardiacmuscle. This tissue-scale model considers the my-ocardium to be composed of continuous intracellularand extracellular compartments, coupled via a contin-uous cellular membrane. The bidomain equations canbe derived from a model that accounts for the tis-sue microarchitecture , but the bidomain model isa fundamentally homogenized description of excitationpropagation that neglects the details of this microarchi-tecture. Instead, the bidomain equations describe thedynamics of a local average of the voltages in the intra-cellular and extracellular compartments over a controlvolume. One of the assumptions required by the homog-enization procedure is that the control volume is largecompared to the scale of the cellular microarchitecturebut small compared to any other important spatial scaleof the system, such as the width of the action potentialwavefront. Although the validity of this model has beenquestioned, for example by Bueno-Orovio et al. , thisapproach has been extremely successful, and at present,most simulations of cardiac electrophysiology use suchmodels. For a detailed review of the bidomain model andother models of electrophysiology, we refer to Griffith andPeskin . Here, we assume that the homogenization as-sumptions hold, and we derive the hyperbolic bidomain V e ( x ) R e ∆ x I e ( x ) L e ∂I e ∂t ∆ x V e ( x + ∆ x ) Extracellular compartment I t ∆ x I ion ∆ xC m ∆ x I t ∆ x I ion ∆ xC m ∆ xV i ( x ) R i ∆ x I i ( x ) L i ∂I i ∂t ∆ x V i ( x + ∆ x ) Intracellular compartmentMembrane
Figure 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). model phenomenologically, starting from charge conser-vation in quasistatic conditions.In the hyperbolic bidomain model, as in the stan-dard bidomain model, the extracellular and intracellu-lar compartments are characterized by the anisotropicconductivity tensors, respectively σ e and σ i . Intro-ducing a local orthonormal basis, { f , s , n } , we as-sume the conductivity tensors can be written as σ j = σ f j f ⊗ f + σ s j s ⊗ s + σ n j n ⊗ n , for j = e , i . Typicalvalues of the conductivity coefficients range from − mS/cm in the cross-fiber direction to 1 mS/cm in thefiber direction (for example, Colli Franzone et al. use σ f = 1 . mS/cm, σ s = 0 . mS/cm, and σ n = 0 . mS/cm). We model the current density fluxes using aCattaneo-type equation, so that τ e ∂ J e ∂t + J e = − σ e ∇ V e , (1) τ i ∂ J i ∂t + J i = − σ i ∇ V i , (2)where V i ( V e ), J i ( J e ), and τ i ( τ e ) are the intracellu-lar (extracellular) potential, the intracellular (extracellu-lar) flux, and the intracellular (extracellular) relaxationtime, respectively. As shown by Jou et al. , relations(1) and (2) can be also derived from the simplest modelof ionic conduction in a dilute system. The derivation ofthese evolution laws for the fluxes from the generalizedGibbs equation shows that they are consistent withthe theory of extended irreversible thermodynamics. Al-ternatively, equations (1) and (2) can be interpreted asarising from a circuit model that includes inductances,such as the one depicted in Fig. 1. The derivation ofthe cable equation for the circuit model shown in Fig. 1is performed in Appendix A. In particular, we show inequation (A15) the relationship between inductance andthe relaxation time.To derive the higher-dimensional model equations, webegin by taking the divergence of the fluxes (1) and (2), so that τ e ∂∂t ∇ · J e + ∇ · J e = −∇ · σ e ∇ V e , (3) τ i ∂∂t ∇ · J i + ∇ · J i = −∇ · σ i ∇ V i . (4)As in similar derivations of the bidomain model, we im-pose a quasistatic form of charge conservation , yielding ∇ · ( J i + J e ) = 0 . (5)Additionally, the current leaving each compartmentneeds to enter the other, so that − ∇ · J i = I t = ∇ · J e , (6)where I t = χ (cid:18) C m ∂V∂t + I ion (cid:19) is the usual transmem-brane current density, with χ the membrane area per unitvolume of tissue, C m the membrane capacitance, and I ion the transmembrane ionic current. Using equation (6) in(3) and (4), we obtain the system of equations, τ e ∂I t ∂t + I t = −∇ · σ e ∇ V e , (7) − τ i ∂I t ∂t − I t = −∇ · σ i ∇ V i . (8)Defining the transmembrane potential V = V i − V e toeliminate V i from the equations yields τ e ∂I t ∂t + I t = −∇ · σ e ∇ V e , (9) τ i ∂I t ∂t + I t = ∇ · σ i ∇ V + ∇ · σ i ∇ V e . (10)Expanding the transmembrane currents, we finally ob-tain the hyperbolic bidomain model, τ e C m ∂ V∂t + C m ∂V∂t + ∇ · D e ∇ V e = − I ion − τ e ∂I ion ∂t , (11) τ i C m ∂ V∂t + C m ∂V∂t −∇ · D i ∇ V − ∇ · D i ∇ V e = − I ion − τ i ∂I ion ∂t , (12)where we set D i = σ i /χ and D e = σ e /χ . Alternatively,the model equations can be written as τ i C m ∂ V∂t + C m ∂V∂t −∇ · D i ∇ V − ∇ · D i ∇ V e = − I ion − τ i ∂I ion ∂t , (13) ( τ e − τ i ) C m ∂ V∂t + ∇ · D i ∇ V + ∇ · ( D e + D i ) ∇ V e =( τ i − τ e ) ∂I ion ∂t . (14)Notice that if τ i = τ e = τ , then these equations re-duce to a hyperbolic-elliptic system that is similar to theparabolic-elliptic form of the standard bidomain model, τ C m ∂ V∂t + C m ∂V∂t −∇ · D i ∇ V − ∇ · D i ∇ V e = − I ion − τ ∂I ion ∂t , (15) ∇ · D i ∇ V + ∇ · ( D e + D i ) ∇ V e = 0 . (16)If we further take τ = 0 , we retrieve the usual bidomainmodel in its parabolic-elliptic form, C m ∂V∂t − ∇ · D i ∇ V − ∇ · D i ∇ V e = − I ion , (17) ∇ · D i ∇ V + ∇ · ( D e + D i ) ∇ V e = 0 . (18) II. REDUCTION TO THE HYPERBOLIC MONODOMAIN
The hyperbolic bidomain model can be simplified byassuming the extracellular and intracellular compart-ments have the same anisotropy ratios, so that D = D e = λ D i . If we make this assumption, we obtain from(14) −∇ · D ∇ V e = 1 λ + 1 ∇ · D ∇ V + λλ + 1 ( τ e − τ i ) C m ∂ V∂t − λλ + 1 ( τ i − τ e ) ∂I ion ∂t , (19)Substituting in (13), we find (cid:20) τ i + λλ + 1 ( τ e − τ i ) (cid:21) C m ∂ V∂t + C m ∂V∂t − λλ + 1 ∇ · D ∇ V = − I ion − (cid:20) τ i + λλ + 1 ( τ e − τ i ) (cid:21) ∂I ion ∂t . (20)Defining τ = τ i + λ ( τ e − τ i ) / ( λ + 1) , we obtain the hy-perbolic monodomain model, τ C m ∂ V∂t + C m ∂V∂t − ∇ · D ∇ V = − I ion − τ ∂I ion ∂t , (21) C m σ f σ s = σ n χ k b µ µ ε where here we have absorbed the term λ/ ( λ + 1) into D . The relaxation time τ of the monodomain model isalways positive, and it is zero only if both τ i and τ e arezero. 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 = ∂V∂t , we transformthe hyperbolic monodomain equations into the first-ordersystem, ∂V∂t = Q, (22) τ C m ∂Q∂t + C m Q − ∇ · D ∇ V = − I ion − τ ∂I ion ∂t . (23)Equations (22)–(23) are usually supplemented with insu-lation boundary conditions, such that ∇ V · N = 0 (24)on the boundary of the domain, where the vector N isthe normal to the boundary.Following Stan et al. , we say that a solution of thehyperbolic monodomain equations has a finite propaga-tion speed if, given compactly supported initial condi-tions for V at time t = 0 , V ( · , t ) is also compactly sup-ported for any t > . The compact support is taken withrespect to the resting potential V . By contrast, a so-lution has infinite propagation speed if the initial dataare compactly supported, but for any t > and any R > , the set M R,t = { x : || x || ≥ R, V ( x , t ) > V } haspositive measure. For any relaxation time τ (cid:54) = 0 , (21)is hyperbolic, and it thereby has the property of finitepropagation speed. Setting τ = 0 , the equations becomeparabolic and have solutions with infinite propagationspeeds. Specifically, in the standard parabolic model, lo-cal perturbations in V , even those that do not generatea propagating front, will travel at infinite speed. III. IONIC MODELS
Transmembrane ionic fluxes through ion channels,pumps, and exchangers are responsible for the cardiac ac-tion potential. The action potential is initiated by a fastinward sodium current that depolarizes the cellular mem-brane. After depolarization phase, slow inward currents(primarily calcium currents) and slow outward currents(primarily potassium currents) approximately balanceeach other, prolonging the action potential and creatinga plateau phase. Ultimately, the slow outward currentsbring the transmembrane potential difference back to itsresting value of approximately − mV. The bidomain C m σ f σ s = σ n χ τ + v τ − v τ − v τ + w τ − w τ d τ τ r τ si k V sic V c V v parameter set 3 1 0.1 0.0125 1 3.33 19.6 1250 870 41 0.25 12.5 33.33 29.0 10 0.85 0.13 0.04parameter set 4 1 0.1 0.0125 1 3.33 15.6 5 350 80 0.407 9.0 34.0 26.5 15 0.45 0.15 0.04parameter set 5 1 0.1 0.0125 1 3.33 12.0 2 1000 100 0.362 5.0 33.33 29.0 15 0.7 0.13 0.04parameter set 6 1 0.1 0.0125 1 3.33 9.0 8 250 60 0.395 9.0 33.33 29.0 15 0.5 0.13 0.04Table II. Sets of parameters for the Fenton-Karma ionic model (27)–(31) as stated by Fenton et al. . and monodomain models must be completed by specify-ing the form of I ion , which accounts for these transmem-brane currents. The states of the transmembrane ionchannels are described by a collection of variables, w , as-sociated with the ionic model, so that I ion = I ion ( V, w ) .Typically, the state variable dynamics are determined bya spatially decoupled system of nonlinear ordinary differ-ential equations, ∂ w ∂t = g ( V, w ) , (25)where the form of g depends on the details of the partic-ular ionic model. We consider the simplified piecewise-linear model of McKean , the two-variable model ofAliev and Panfilov , the three-variable model of Fentonand Karma , and the biophysically detailed models often Tusscher and Panfilov for the ventricles (20 vari-ables) and of Grandi et al. for the atria (57 variables).In the two-variable model of Aliev and Panfilov , w = { r } , with ∂r∂t = (cid:18) ε + µ rµ + V (cid:19) ( − r − kV ( V − b − , (26)and the ionic current is I ion = kV ( V − α ) ( V −
1) + rV .The parameters used in this case are shown in Table I. Inthe three-variable model of Fenton and Karma , w = { v, w } , with ∂v∂t = 1 τ − v ( V ) (1 − p ) (1 − v ) − τ + v pv, (27) ∂w∂t = 1 τ − w (1 − p ) (1 − w ) − τ + w pw, (28)where τ − v ( V ) = (1 − q ) τ − v + qτ − v , p = H ( V − V c ) , q = H ( V − V v ) , and H ( · ) is the Heaviside function. Thetotal ionic current is the sum of three currents, I ion = − I fi ( V, v ) − I so ( V ) − I si ( V, w ) , with I fi ( V, v ) = − τ d vp ( V − V c ) (1 − V ) , (29) I so ( V ) = 1 τ V (1 − p ) + 1 τ r p, (30) I si ( V, w ) = − τ si w (cid:0) (cid:0) k (cid:0) V − V sic (cid:1)(cid:1)(cid:1) . (31)The sets of parameters used in this paper for this modelare shown in Table II. We refer to the original papers for the statements of the equations and parameters of thebiophysically detailed ionic models. IV. NUMERICAL RESULTS
The hyperbolic monodomain and bidomain models arediscretized using a low-order finite element scheme de-scribed in Appendix C that is implemented in the open-source parallel C++ code BeatIt (available at http://github.com/rossisimone/beatit ), which is based onthe libMesh finite element library and relies on lin-ear solvers provided by PETSc . All the code usedfor the following tests is contained in the online reposi-tory, and all tests can be replicated directly from thosecodes. The only exception is the atrial fibrillation test,which uses a patient-specific mesh that is not containedin the repository. One- and two-dimensional simulationswere run using a Linux workstation with two Intel XeonE5-2650 v3 processors (up to 40 threads) and 32 GB ofmemory. Three-dimensional simulations were run on theKillDevil Linux cluster at the University of North Car-olina at Chapel Hill. We used Matlab to visualize theone-dimensional results and Paraview for the two- andthree-dimensional simulations. A. Comparison with an exact solution
We start by considering a simple piecewise-linearbistable model for the ionic currents, with I ion ( V ) = kV − k [ V H ( V − V ) + V (1 − H ( V − V ))] , (32)where V is the resting potential, V is the threshold po-tential, V is the depolarization potential, and H ( · ) isthe Heaviside function. This model (32) reduces to thepiecewise-linear model of McKean after a simple di-mensional analysis, after which the nondimensional ioniccurrents take the form ˆ I ion (cid:16) ˆ V (cid:17) = ˆ V − H (cid:16) ˆ V − α (cid:17) . (33)The nondimensional potential ˆ V is related to the dimen-sional potential by V = ( V − V ) ˆ V + V . In this model, α = ( V − V ) / ( V − V ) describes the excitability of thetissue. As we show in Appendix B, using model (33),it is possible to find the analytic solution of a propagat-ing front for the nondimensional hyperbolic monodomainproblem. In particular, we find that the front propaga-tion speed in an unbounded domain is c = (1 − α ) (cid:113) µ + ( α − α ) ( µ − < (cid:114) µ = c s . (34) [ms] nond i m e n s i on a l c ondu c ti on v e l o c it y exact = 0.1numerical = 0.1exact = 0.2numerical = 0.2exact = 0.3numerical = 0.3c s x = 0 = 0.1 = 0.3 x = 0 = 0.1 = 0.3 Figure 2. Conduction velocities for the piecewise-linear reaction model (33). Left) comparison between the numerical and theexact wave speed (34). The black curve represent the characteristic propagation speed of the system, i.e. the limit speed atwhich 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 fordifferent relaxation times at t = 10 . All simulations were run using a mesh size h = 0 . and a time step size ∆ t = 0 . . The speed c s = (cid:112) /µ = (cid:112) C m /τ k represents the maxi-mum speed at which a perturbation can travel in the sys-tem. When the relaxation time τ goes to zero, the localperturbations can travel at infinite speed. The parame-ter µ = τ k/C m is a nondimensional number representingthe ratio between the relaxation time and the character-istic time of the reactions and characterizing the effectsof the inductances in the system (see Appendix B). Thecorresponding dimensional speed is v = (cid:115) σkχC m (1 − α ) (cid:113) µ + ( α − α ) ( µ − < (cid:114) σχτ C m . (35)As a verification test, we consider the spatial interval Ω = [0 , . An initial stimulus of amplitude is ap-plied at x ∈ [24 . , . for the interval t ∈ [0 . , . .The system of equations is solved using a mesh size h = 0 . and a time step size ∆ t = 0 . . Weshow the nondimensional solutions ˆ V and ˆ Q at t = 10 in Fig. 2(center and right). We register the activationtime at a particular spatial location whenever ˆ V crossesa threshold of . . The conduction velocities are mea-sured by picking the distance between two points x and x ∈ Ω and dividing it by the time interval between theactivation times at these two locations. Because (34)is obtained by assuming that Ω is the entire real line,we need to measure conduction velocities far from theboundaries. For this reason, we choose x = 30 and x = 32 . The comparison between the exact conductionvelocities defined by equation (34) and those obtained bythe simulations is shown in Fig. 2(left).We also examine the effect of the relaxation time onthe conduction velocities for three values of the excitabil-ity parameter, α ∈ { . , . , . } . With this simple ionicmodel, the larger the relaxation time, the slower the wavespeed. The limit speed c s at which local perturbationscan travel is also shown in Fig. 2(left). Thus, in thispiecewise-linear model, the effect of inductances is to slowdown the propagation of the front. By contrast, the fol- lowing tests will show that with fully nonlinear models,inductances can also enhance propagation.Using the analytic solution derived in Appendix B, wealso perform a space-time convergence study. We con-sider the spatial interval Ω = [ − , and set the initialconditions according to (B9) and (B10), assuming thefront at time t = 0 is at x = 0 and α = 0 . . The system isdiscretized in time using implicit-explicit (IMEX) Runge-Kutta (RK) time integrators . In particular, we usethe forward/backward Euler schemes for the first-ordertime integrator and the Heun/Crank-Nicholson schemesfor the second-order time integrators. For more detailson the time integrator, see Appendix C. To ensure thevalidity of the solutions (B9) and (B10) in the consid-ered 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 proportionalto the mesh size h , with a value of ∆ t = 0 . in the coars-est cases. The time step size was chosen to guarantee alarge enough number of time iterations and while satisfy-ing the CFL condition (see Appendix C). Fig. 3 shows theerrors for the potential V and its time derivative Q usingthe first-order (Fig. 3 left) and second-order (Fig. 3 cen-ter) time stepping schemes. Because the ionic currentsand their derivative in this case are not smooth func-tions, we do not expect to obtain full second-order con-vergence. On the other hand, whereas the convergencerates for the parabolic monodomain model are alwaysfirst order in both V and Q , the hyperbolic monodomainmodel with the second-order time stepping scheme con-verges quadratically in V and linearly only in Q . Thesimulations were run with and without regularization ofthe Heaviside and Dirac- δ functions. The errors reportedin Fig. 3 correspond to the more accurate simulations.We show in Fig. 3(right) that the slow convergence re-sults from the non-smoothness of the ionic currents: wecompare the convergence of the piecewise-linear reactionmodel with the monodomain model with the cubic reac- -3 -2 -1 h -4 -2 || V - V e x ac t || L = 0.6 ms = 0.4 ms = 0.2 ms = 0 msh -3 -2 -1 h -6 -4 -2 || V - V e x ac t || L = 0.6 ms = 0.4 ms = 0.2 ms = 0 mshh -3 -2 -1 h -5 || V - V e x ac t || L piecewise linear = 0 mscubic = 0 mshh -3 -2 -1 h -3 -2 -1 || Q - Q e x ac t || L = 0.6 ms = 0.4 ms = 0.2 ms = 0 msh -3 -2 -1 h -4 -2 || Q - Q e x ac t || L = 0.6 ms = 0.4 ms = 0.2 ms = 0 msh -3 -2 -1 h -10 -5 || Q - Q e x ac t || L piecewise linear = 0 mscubic = 0 mshh Figure 3. Convergence study for the piecewise-linear reaction model (33) for the variables V and Q . Left) comparison of the rateof 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 usinga second-order IMEX-RK time integrator. Right) comparison of the rate of convergence between the piecewise-linear reactionmodel and a cubic reaction model for the parabolic monodomain model using a second-order IMEX-RK time integrator. tion term ˆ I ion ( ˆ V ) = ˆ V ( ˆ V − V − α ) . The exact solution for the cubic reaction in the parabolicmonodomain has been determined previously . Weare not aware of the existence of an analytical solutionfor the hyperbolic monodomain with a cubic reactionterm. We show in Fig. 3 that the cubic model con-verges quadratically (using a second-order time steppingscheme) whereas the piecewise-linear model convergeslinearly. To conclude, at least in these tests, the dis-cretization errors of the parabolic mondomain model arelarger than the discretization errors in the hyperbolicmonodomain model, indicating a better accuracy for thehyperbolic model.
B. Effect of the relaxation time on the conduction velocityfor different ionic models
We now analyze how the inductance terms influencethe conduction velocities of cardiac action potentials. Weconsider four ionic models: the two-variable model ofAliev and Panfilov , the three-variable model of Fentonand Karma , the ventricular M-cell model of ten Tuss-cher and Panfilov (20 variables), and the atrial modelof Grandi et al. (57 variables). For the two-variableAliev-Panfilov model, we use the parameters of Nash and Panfilov , as reported in Table I. We consider two val-ues for the excitability parameter, α = 0 . and α = 0 . .For the three-variable Fenton-Karma model, we considerparameters sets 3, 4, 5, and 6 given by Fenton et al. and reported in Table II. For the biophysically detailedionic model of Ten tusscher et al. and Grandi et al., oursimulations used the C++ implementations of the mod-els provided by the authors, using C m = 1 µ F / cm forthe membrane capacitance in both models.For all cases, the domain is the interval Ω = [0 , cm.The simulations use a mesh size h = 31 . µ m and a timestep size ∆ t = 0 . ms. An initial stimulus is appliedat the center of the domain, x ∈ [2 . , . cm, duringthe time interval t ∈ [0 . , . ms. The conductionvelocities are measured, as in the previous test case, bydividing the distance between two points x and x by thetime interval between the activation times at these twoselected locations. The computed conduction velocitiesof different ionic models are shown in Fig. 4 for relaxationtimes in the range [0 , ms. Contrary to the results forthe piecewise-linear ionic model, in this case, small andmoderate values of the relaxation time act to enhancepropagation.We also compare the computational time required bythe parabolic and hyperbolic monodomain models in twoand three spatial dimensions. Table III shows the numberof iterations taken by the linear solvers (conjugate gra-dient preconditioned by successive over-relaxation) alongwith the relative computational times of the hyperbolic [ms] nond i m e n s i on a l c ondu c ti on v e l o c it y = 0.1 = 0.2 [ms] c ondu c ti on v e l o c it y [ c m / s ] param set 3param set 4param set 5param set 6 [ms] c ondu c ti on v e l o c it y [ c m / s ] Grandi model '11Ten tusscher model '06
Figure 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 thecenter of a one-dimensional domain and the activation times are recorded. For all models we used a time step size ∆ t = 0 . ms and a mesh size h = 31 . µ m over a domain of length cm. Unlike the linear case, the nonlinearities of the time derivativeof the ionic currents increase the wavefront speed for small and moderate values of the relaxation time τ . monodomain model (with respect to the parabolic model)for the reaction step and the diffusion step . In the reac-tion step, we solve the ionic model, and we also evaluatethe ionic currents and their time derivatives. For this rea-son, we expect this step to take longer in the hyperbolicmonodomain model. In fact, with simple ionic models,the amount of work required to evaluate the ionic cur-rents and their time derivatives is almost doubled. This isnatural because, for such simple models, the time deriva-tive of the ionic currents do not use any quantities al-ready computed in the evaluation of the ionic currents.By contrast, for biophysically detailed models, most ofthe computation is needed for the evaluation of the ioniccurrents, and many quantities can be reused for the eval-uation of their time derivatives. This is reflected in TableIII by a small increase (less than 7%) in the computa-tional time of the hyperbolic model for the Ten Tusscher’06 model. Unexpectedly, the solution of the linear sys-tem take fewer iterations in the hyperbolic monodomainmodel. This is reflected by a speedup of about 20% ofthe computational time in the diffusion step, where weassembled the right hand side, solve the linear system,and update the variables Q and V . C. Conduction velocity anisotropy
This section analyzes whether inductances affect theanisotropy in conduction velocity. We define theanisotropy ratio as the ratio between the conductionvelocities in the transverse and longitudinal directions,and to obtain a simple estimate for this ratio, we con-sider plane-wave solutions propagating in these direc-tions. We use the Fenton-Karma model with parame-ter set 3 and consider four values of the conductivity: σ = 0 . mS/mm, σ = 0 . mS/mm, σ = 0 . mS/mm, and σ = 0 . mS/mm. Evaluating the wavespeeds v , v , v , and v corresponding to the conduc-tivities σ , σ , σ , and σ , respectively, we show how the anisotropy ratio changes with respect to the relaxationtime for several conductivity ratios. In fact, as shownin Table II, σ and σ are the typical longitudinal andtransverse conductivities for this model. Considering theinterval Ω = [0 , cm, we use a mesh size h = 25 µ m anda time step size ∆ t = 0 . ms. Fig. 5(left) shows howthe velocities at different σ compare to each other. Tomake the comparison more clear, we normalize the resultswith respect to the conduction velocities of the parabolicmonodomain model, i.e., corresponding to τ = 0 ms. Theresults shown in Fig. 5(left) should be understood as thepercentage difference of the wave speed with respect tothe parabolic case. It is clear that the curves are similarfor the conductivities considered. On the other hand, therelaxation time at which the conduction velocity of thehyperbolic monodomain model matches the conductionvelocity of the parabolic monodomain model decreaseswith smaller conductivities. Fig. 5(right) shows the re-laxation time needed to maintain the same conductionvelocity as in the parabolic monodomain model. In par-ticular, for σ , we find that the relaxation time needed inthe hyperbolic monodomain model to yield the same ve-locities as in the parabolic monodomain models is about τ = 0 . ms. The differences in the velocities computedwith τ = 0 . ms and τ = 0 . ms are smaller than 2%.This difference is typically smaller than the numerical er-ror in the simulations. For this reason, in the followingtests, our comparison will only use the value τ = 0 . ms. Fig. 5(center) shows how the anisotropy ratio isinfluenced by the relaxation time. Although the ratiois not constant, the variation over the relaxation timesconsidered never exceeds 5%. For the conductivity ra-tios 8:1 ( σ : σ ), 4:1 ( σ : σ ), and 2:1 ( σ : σ ), we find thatthe ratios between the conduction velocities in the trans-verse and in the longitudinal directions are approximately1/ √ ( v / v ), 1/ ( v / v ), and 1/ √ ( v / v ). This be-havior is expected because the conduction velocity de-pends on the square root of the conductivity. We con-clude that the effect of the inductances on the anisotropy
2D Max CG Iterations Reaction Step Cost Diffusion Step Cost Overall100 × τ = 0 ms τ = 0 .4 ms τ = 0 . ms /τ = 0 ms τ = 0 . ms /τ = 0 ms τ = 0 . ms /τ = 0 msAP 8 3 -54% 23% -3%FK 22 15 -62% 22% -2%TP06 11 5 -7% 15% -3%3D Max CG Iterations Reaction Step Cost Diffusion Step Cost Overall10 × × τ = 0 ms τ = 0 .4 ms τ = 0 . ms /τ = 0 τ = 0 . ms /τ = 0 ms τ = 0 . ms /τ = 0 msAP 43 16 -73% 48% 30%FK 43 16 -49% 47% 25%TP06 43 16 -6% 46% 15%Table III. Comparison of the compuational cost for the parabolic ( τ = 0 ms) and hyperbolic ( τ = 0 . ms) monodomain modelsfor two- and three-dimensional problems in a serial computation. Negative values indicate relative slowdown of the hyperbolicmodel compared to the parabolic model, whereas positive values represent relative speedup. In the reaction step , we measurethe time needed to solve the ionic model and to evaluate the ionic currents and their time derivatives nodewise for 1000 timesteps. In the diffusion step , we measure the time needed to form the right hand side, solve the linear system and update thesolutions, for 1000 time steps. The hyperbolic monodomain model can be solved using fewer linear solver iterations. This isreflected 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. Forthis reason, the reaction step is actually slower than in the parabolic monodomain model. Notice, however, the difference inthe 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. [ms] no r m a li ze d w a v e s p ee d = 0.1 = 0.05 = 0.025 = 0.0125 [ms] a n i s o t r opy r a ti o / / / [S/mm] r e l a x a ti on ti m e [ m s ] Figure 5. Variation of the conduction velocity with respect to the conducitivity coefficients using the Fenton-Karma model andparameter set 3 in Table II. Left) wave speed for different conductivities, normalized with respect to the parabolic conductionvelocity. 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 aredifferent for different conductivities, the error in the conduction velocity we commit considering the fixed value τ = 0 . ms issmaller than 2% and therefore typically smaller than the spatial error. ratio is negligible. D. Discretization error: a simple two-dimensional benchmark
It is important to be aware of the limitations of thenumerical method used in the simulations. For example,spiral break up can occur easily if the wavefront is not ac-curately captured. In fact, numerical error can introducespatial inhomogeneities that lead both to spiral wave for-mation and also to spiral break up. This numerical arti-fact disappears under grid refinement, and on sufficientlyfine computational meshes, spiral wave break up resultsonly 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 , × [0 , cm, and we use the Fenton-Karmamodel with the parameter set 3. Reentry is induced byan S1-S2 protocol. First, an external stimulus of unitarymagnitude is applied in the left bottom corner, Ω stim = { x ∈ Ω : (cid:107) x (cid:107) ≤ } at t = 0 ms. A second stimulus withthe same amplitude and in the same region is appliedafter 300 ms. We consider two cases: 1) the fiber field isaligned with some of the mesh edges; and 2) the fiber fieldis not aligned with any edge in the mesh. Figs. 6(top) and7(top) demonstrate that this can be easily achieved byfixing the fiber field and rotating the mesh. Those figuresshow the fiber field in green on top of the computationalgrid. The red region at the bottom left of the domain isthe stimulus region, Ω stim . In the first test (Fig. 6), the0 τ = 0 ms h ≈ µ m h ≈ µ m τ = 0 . ms h ≈ µ m h ≈ µ mFigure 6. Effect of fiber direction and mesh orientation onpropagation. Top) fibers field (green), mesh orientation, andstimulus region (red); Left) Monodomain model; Right) Hy-perbolic monodomain model. Domain size is 12x12 cm andthe solution is evaluated at time t = 450 ms. fibers are set to be orthogonal to the propagating front.Using 512 elements per side, which is equivalent to amesh size h of approximately 234 µ m (163 µ m in thedirection of wave propagation), the solutions obtainedon the two meshes differ substantially. Large differencescan 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). Itis clear from Fig. 6 that whenever the fiber field is notaligned with the mesh, the conduction velocity is largelyoverestimated. The second test is similar to the first one,but the fibers are now rotated by 90 ◦ ; see Fig. 7. It isclear that the error on the conduction velocity in the fiberdirection is much smaller than the error in the transversedirection. In fact, the solutions for h ≈ µ m and h ≈ µ m show the greatest differences for transversepropagation.Although similar results have already been reported inprior work , it is nonetheless generally believed thata mesh size of the order of 200 µ m is usually sufficient forcardiac computational electrophysiology . This beliefis supported by numerical simulations using a mesh sizeof 200 µ m that show that the error on the wave speed inthe fiber direction is less than 5%. On the other hand,the most interesting dynamics of the propagating frontoften will occur perpendicular to the fiber direction. Forexample, during the normal electrical activation of theventricles, the signal spreads from the endocardium tothe epicardium traveling across the ventricular wall, per-pendicular to the fibers. The reduced conductivity in thetransversal direction, usually taken to be about 8 timessmaller, requires a finer resolution of the grid. As notedby Quarteroni et al. , the required grid resolution to τ = 0 ms h ≈ µ m h ≈ µ m τ = 0 . ms h ≈ µ m h ≈ µ mFigure 7. Effect of fiber direction and mesh orientation onpropagation. Top) fibers field (green), mesh orientation, andstimulus region (red); Left) Monodomain model; Right) Hy-perbolic monodomain model. Domain size is 12x12 cm andthe solution is evaluated at time t = 360 ms. capture the transverse conduction velocity with an errorsmaller than 5% is about 25 µ m. These results stronglyindicate that to capture correctly the wave speeds, themesh discretization needs to be about the size of a singlecardiomyocyte, as also discussed by Hand and Griffith . E. Effect of the relaxation time on spiral break up
Our next tests explore the effect of relaxation time onspiral wave break up. We consider a square slab of tissue,
Ω = [0 , × [0 , cm. An initial stimulus is applied inthe region y < . cm for 1 ms at t = 0 ms to generatea wave propagating in the y -axis. A second stimulus isthen applied in the region { x < cm ∧ y < cm } for 1ms at time t = 320 ms to initiate a spiral wave. Spiralbreak up is easily obtained using the parameter set 3 inTable II, because of the steep action potential duration(APD) restitution curve. In this case, the back of thewave forms scallops, and when the turning spiral tries toinvade these regions, it encounters refractory tissue andbreaks.Fig. 8 shows the formation of the spiral wave in thecase of isotropy. With τ = 0 . ms and the Fenton-Karmamodel with parameter set 3, the conduction velocity ofthe hyperbolic monodomain model is very close to theconduction velocity of the parabolic monodomain model.In fact, the evolution of the spiral waves in the two casesis very similar.Fig. 9 compares spiral break up obtained using themonodomain model and the hyperbolic monodomainmodel in the anisotropic case, with fibers aligned withthe y -axis. The relaxation time for the hyperbolic mon-1 τ = 0 ms τ = 0 . ms V Q V Q
Time: 400 msTime: 500 msTime: 600 msTime: 700 msTime: 800 msFigure 8. Comparison of the evolution of a spiral wave forthe isotropic monodomain model (left) and the isotropic hy-perbolic monodomain model (right) using the Fenton-Karmamodel with parameter set 3. The wavefronts (red) and thetails (light blue) represented by the variable Q are shownnext to the transmembrane potential V . The evolution ofthe spiral wave is similar in the two cases. The domain is [0 , × [0 , cm, the mesh size is h ≈ µ m, and thetime step size ∆ t = 0 . ms. The fibers are aligned with the y -axis. odomain model is chosen to be τ = 0 . ms. As explainedin Sec. IV C, although the conduction velocities in thetransverse direction are different in the hyperbolic andparabolic models for τ = 0 . ms, the difference is smallerthan 2%. This difference is smaller than the spatial errorand, for this reason, we consider the two cases to yieldessentially the same anisotropy ratio. It is clearly shownin Fig. 9 that in the initial phase (up to t = 600 ms), thespirals are almost identical. After the first rotation, how-ever, the spiral wave starts to break, and the relaxationtime of the hyperbolic monodomain model shows quan-titative differences in the form of the break up. Fig. 9shows the dynamics of Q = ∂ t V . This variable can beused to define the fronts (red) and the tails (light blue)of the waves. τ = 0 ms τ = 0 . ms V Q V Q
Time: 600 msTime: 800 msTime: 1000 msTime: 1250 msTime: 1500 msFigure 9. Comparison of the evolution of a spiral break upfor the monodomain model (left) and the hyperbolic mon-odomain model (right) using the Fenton-Karma model withparameter set 3. The wavefronts (red) and the tails (lightblue) represented by the variable Q are shown next to thetransmembrane potential V . Although the initiation sequenceis the same and the spiral wave at t = 600 ms are similar,the evolution of the break up is different. The domain is [0 , × [0 , cm, the mesh size is h ≈ µ m in the fiberdirection, and h ≈ µ m in the transverse direction, and thetime step is ∆ t = 0 . ms. The fibers are aligned with the y -axis. F. Atrial fibrillation
As a more complex application of the model, we usethe hyperbolic monodomain model to simulate atrial fib-rillation. Although a rigorous study of atrial fibrillationrequires the use of a bidomain model, similar to the oneproposed in equations (11) and (12), the results given bythe monodomain model can represent a reasonable ap-proximation to the dynamics of the full bidomain model.For instance, Potse et al. provide a comparison of thedynamics of the bidomain and monodomain models.The anatomical geometry used in these simulationswas based on the human heart model constructed2 Rule-based fiber field τ = 0 [ ms ] τ = 0 . ms ] Figure 10. Left) Reconstructed fiber field of left atria based on the data in . Right) Activation times for the monodomain andfor the hyperbolic monodomain model with τ = 0 . ms using the Fenton-Karma model with parameter set 3. τ = 0 [ ms ] Time: 190 ms 420 ms 980 ms 1980 ms τ = 0 . ms ] Time: 190 ms 420 ms 980 ms 1980 msFigure 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 . ms. Using the hyperbolicmonodomain model spiral waves form after second stimulus. The parabolic monodomain model requires 3 stimuli before spiralwaves 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 . by Segars et al. using a 4D extended cardiac-torso(XCAT) phantom. From their data, we extracted andreconstructed a geometrical representation of the leftatrium using SOLIDWORKS. We then used Trelis togenerate a simplex mesh of the left atrium consisting ofapproximately 3.5 million elements. Our simulations usethe Fenton-Karma model with parameter set 3.A major challenge in modeling the atria is the defi-nition of the fiber field. In work spanning the past 100years, several studies have tried to characterize the fiberarchitecture of the atria , but the structure of themuscle in the atria is so complex (as can be seen from thediagrams by Thomas ) that developing a set of math-ematical rules to reproduce a realistic fiber field is chal-lenging. Our strategy for modeling the fiber field is touse the gradient of a harmonic function that satisfies pre-scribed boundary conditions. We divided the left atriuminto several regions, and solved three Poisson problems:one for the left atrial appendage; one for the Bachmann’sbundle; and one for the remaining parts of the left atrium.In particular, for the left atrial appendage, we used amethod similar to the one proposed by Rossi et al. In the other regions, similar to work by Patelli et al. andKrueger et al. , we used the normalized gradient of thesolution of the Poisson problem to reconstruct a fiberfield similar to the one shown in the studies by Ho etal. . The approximate fiber field of the left atrium re-constructed using this strategy is shown in Fig. 10, wherethe colors represent the magnitude of the x -componentof the fibers and are used only to highlight changes indirection.To initiate atrial fibrillation, we use an S1-S2 stimu-lation protocol that is applied at the junction with theinter-atrial band, as described by Colman . The S1stimulus is applied with a cycle length of 350 ms fol-lowed by a short S2 stimulus with a cycle length of 160ms. Fig. 10(right) shows that the initial electrical ac-tivation times are similar for the monodomain and hy-perbolic monodomain models. Because this set of pa-rameters gives a velocity of about 45 cm/s, roughly 35%slower than the expected conduction velocity , it takeslonger for the atrium to be fully activated. When thesecond stimulus is applied, spiral waves form in the hy-perbolic monodomain system but not in the parabolic3 τ e = 0 ms, τ i = 0 ms τ e = 0 . ms, τ i = 0 ms τ e = 0 . ms, τ i = 0 ms τ e = 0 ms, τ i = 0 . ms τ e = 0 . ms, τ i = 0 . ms τ e = 0 . ms, τ i = 0 . ms τ e = 0 ms, τ i = 0 . ms τ e = 0 . ms, τ i = 0 . ms τ e = 0 . ms, τ i = 0 . msFigure 12. Transmembrane potential patterns formed by the hyperbolic bidomain model with the Fenton-Karma model at 2ms after a unipolar cathode stimulus is supplied in the extracellular compartment. The red regions are areas of depolarizationand blue regions are areas of hyperpolarization. The hyperbolic bidomain model shows a behavior qualitatively similar to theconventional bidomain model shown in the study by , but the details of the are influenced by the extracellular and intracellularrelaxation times. monodomain model; see Fig. 11 at 420 ms. Spiral wavesare generated for τ = 0 ms only after the third stimulusis applied. Fig. 11 also shows the fronts (red) and tails(light blue) of the waves, as indicated by Q . G. Virtual electrodes
The reduction of the bidomain model to the mon-odomain model is not valid if the extracellular and intra-cellular compartments have different anisotropy ratios.As shown in experimental studies , when a current issupplied by an extracellular electrode, adjacent areas ofdepolarization and hyperpolarization are formed in thecase of unequal anisotropy ratios. Applying a stimula-tion current I e,stim in the extracellular compartment, theregion of hyperpolarization near a cathode ( I e,stim < )is called a virtual anode, while the region of depolariza-tion near an anode ( I e,stim > ) is called a virtal cathode.The existence of virtual cathodes and anodes is impor-tant to understand the four mechanisms of cardiac stim-ulation and to study the mechanisms of defribillation. Inthis section, we focus only on a unipolar cathodal stim-ulation to investigate the cathode formation mechanismin the hyperbolic bidomain model. We refer the readerinterested in this virtual electrode phenomenon to themore detailed numerical studies of Wikswo and Roth and Colli Franzone et al. .Following Sepulveda et al. , we consider the domain Ω = [ − , cm × [ − . , . cm, and consider this to cor-respond to the epicardial surface. We apply a unipolarcathodal current in the region Ω stim = [ − . , . mm × [ − . , . mm. We use the Fenton-Karma model withparameter set 3, and the corresponding nondimension-alized current stimulus amplitude is − σ fe = 1 . , σ se = σ ne = 1 . and σ fi = 2 . , σ si = σ ni = 0 . , with the fibers aligned with the x -axis.Fig. 12 shows the pattern formed by the transmem-brane potential 2 ms after the extracellular stimulus isapplied. Red regions are areas of depolarization, andblue regions are areas of hyperpolarization. The depo-larization and hyperpolarization regions are affected bythe choice of the extracellular and intracellular relaxationtimes. For this test, we consider the relaxation times τ i and τ e to take the values 0, 0.2, and 0.4 ms. Note thatif τ i = τ e , the hyperbolic bidomain equations (13) and(14) reduce to the hyperbolic-elliptic system of equations(15) and (16). Although this test shows the ability of thehyperbolic bidomain model to reproduce the virtual elec-trode phenomenon, the influence of the relaxation timesin the stimulation remains unclear and necessitates fur-ther investigation.4 V. CONCLUSIONS
Local perturbation in the standard bidomain and mon-odomain models propagate with an infinite speed. Tocorrect this unrealistic feature, we developed a hyperbolicversion of the bidomain model, and in the case that intra-cellular and extracellular conductivity tensors have equalanisotropy ratios, we reduced this model to a hyperbolicmonodomain model. Our derivation relies on a Cattaneo-type model for the fluxes, described by equations (1) and(2). The Cattaneo-type fluxes are equivalent to intro-ducing self-inductance effects, as shown by the schematicdiagram in Fig. 1. As shown in Appendix A, relaxationtimes introduced in the Cattaneo fluxes represent the ra-tio between the inductances and the resistances of theextracellular and intracellular compartments. Althoughthe hyperbolic monodomain model reduces to the classi-cal parabolic model in the case that the relaxation timesof both the intracellular and extracellular compartmentsare zero, the models differ if at least one of these relax-ation times is nonzero.Although hyperbolic bidomain and monodomain mod-els do not appear to have been previously discussed inthe literature, the work of Hurtado et al. does addressthe problem of propagation with infinite speed. Theirapproach uses a nonlinear model for the fluxes based onporous medium assumptions. Because of the nonlinearnature of the fluxes, the porous medium approach to car-diac electrophysiology is difficult to analyze, even for sim-ple linear reactions. Hurtado et al. use the simplifiedionic model proposed by Aliev and Panfilov . It would beinteresting to see the extension to biophysically detailedmodels. Moreover, it is clear that the porous mediumapproach is not incompatible with the theory developedhere. In fact, the two models could be combined using anonlinear Cattaneo-type model for the fluxes.As discussed by King et al. , abnormal conductionvelocities may contribute to arrhythmogenesis. We havestudied how the conduction velocities of the propagat-ing action potential change with respect to the relax-ation time of the Cattaneo-type fluxes. In the case of thepiecewise-linear model of McKean , it is possible to findan analytical expression for the conduction velocity ofthe hyperbolic monodomain model. (Unfortunately, thewave speed of this bistable model has been reported in-correctly in some prior studies .) We have shown thatfor the hyperbolic monodomain model, signals cannottravel faster than the characteristic propagation speedof the medium. Our numerical results match the theo-retical predictions for different values of the excitabilityparameter. In this simple case, conduction velocity is amonotone decreasing function of relaxation time. Addi-tionally, we have shown that, for this linear case, dis-cretization errors are lower in the hyperbolic model.Using simple nonlinear ionic models , however, wefound that the relationship between the relaxation timeand the conduction velocity is not always monotone. Infact, for small and moderate values of the relaxation time, we found that waves propagate faster than in the stan-dard monodomain model. With these ionic models, amaximum conduction velocity is found at an optimalrelaxation time. At small relaxation times, the secondderivative terms are relatively small. Consequently, thisphenomenon likely results from the time derivative ofthe ionic currents. For larger values of the relaxationtimes, inertial effects dominate, and the conduction ve-locity monotonically decreases to zero. This implies thatwe can find a value of the relaxation time for which theaction potential propagates at the same velocity as in thestandard monodomain model. This fact allowed us to di-rectly investigate the influence of inductances in the mon-odomain model. In the biophysically detailed ionic mod-els for ventricular myocytes and atrial myocytes con-sidered here, a similar behavior can be found, althoughonly for very small values of the relaxation time. Forthis reason, we compared the parabolic and hyperbolicmonodomain models in two and three spatial dimensionsusing the Fenton-Karma model .Changing the type of equation may have importantconsequences in terms of the stability of the numericalapproximations. We did not encounter any substantialnumerical difficulties beyond those already present in thestandard parabolic models. On the contrary, we haveshown that the linear system is easier to solve and thesolutions are more accurate when the hyperbolic systemis used. Moreover, we believe that the solution of theparabolic monodomain and bidomain model already havemany of the numerical difficulties usually associated withhyperbolic systems: the propagating action potential inthe parabolic models already have sharp fronts which rep-resent one of the main difficulties in hyperbolic systems.Additionally, the hyperbolic monodomain and bidomainmodels are damped wave equations, where the dampingterm is dominant. In conclusion, the introduction of hy-perbolic terms to the monodomain and bidomain modelsdo not appear to add substantial numerical difficulties.We performed several qualitative comparisons betweenthe parabolic and the hyperbolic monodomain models.In particular, we considered the case of spiral break upresulting from conduction block caused by a steep APDrestitution curve. Using parameter set 3 of the Fenton-Karma model, the back of the wave easily forms a series ofindentations or scallops. When a spiral wave invades thisafter turning, it encounters refractory tissue that causesconduction block, leading to a break up of the wave. Thisbehavior, already observed for the standard monodomainmodel and explained by Fenton et al. , is still presentin the hyperbolic monodomain model. We showed spi-ral break up in a simple two-dimensional test case andalso in an anatomically realistic three-dimensional modelof the left atrium. We reconstructed on the left atria afiber field qualitatively in accordance with the data re-ported by Ho et al. , Krueger et al. , and Pashakhan-loo et al. We applied a fast-pacing S1-S2 stimulationprotocol to the left atrium to induce spiral waves and spi-ral break up. Although the initial activation sequences5were similar, for the hyperbolic monodomain model, spi-ral waves appeared already after the second stimulus wasapplied, whereas for the parabolic monodomain at leastthree stimuli were needed.We also used a simple two-dimensional benchmark toinvestigate how the alignment of the fiber direction withrespect to the elements of the finite element mesh influ-ences wave propagation. Under mesh refinement, spatialdiscretization errors become smaller and smaller, but thecommon assumption that a grid resolution of 200 µ m issufficient to correctly capture the conduction velocity seems to be overly optimistic. In our tests, we analyzedthe propagation of a simple wave using an anisotropicconductivity tensor. Fixing the fiber field and rotatingthe mesh, we were able to show that using a mesh sizeof about 234 µ m (163 µ m in the direction of wave prop-agation), the propagation was greatly influenced by themesh orientation relative to the fiber alignment. Even ona finer mesh with a mesh size of about 117 µ m (82.5 µ min the direction of wave propagation), visible differencesbetween the propagation patterns obtained with the dif-ferent orientations remained. These tests highlight thefact that if the transverse conductivity coefficients aretaken to be about 8 times smaller than the longitudi-nal conductivity, as is typically done in practice ,the mesh size needs to be much smaller than 200 µ m inorder to yield resolved dynamics. Our results are in ac-cordance with the convergence study by Rossi et al. for propagation in the transverse direction, in which theauthors showed that the mesh size in the transverse di-rection needs to be about 25 µ m. Because such meshresolutions require elements only slightly larger than theactual cardiac cells, it is natural to ask whether the useof a continuum model is even appropriate: for a similarcomputational cost, it could be possible to construct adiscrete model of the heart. Although the difficulties inrepresenting correctly the conduction velocities are welldocumented in the literature and have been thor-oughly analyzed by Pezzuto et al. , the focus generallyseems to have been on the conduction in the fiber direc-tion, as is clear from the choice of the three-dimensionalbenchmark test proposed by Niederer et al. Other solu-tion strategies for the monodomain and bidomain mod-els using adaptive mesh refinement and high-orderelements have been proposed, but their application sofar appears to be relatively limited.Finally, we showed that the hyperbolic bidomain modelcan capture the virtual electrode phenomenon observedin the standard model. By stimulating the epicardial sur-face using a unipolar cathodal current, we have shownhow the pattern of the transmembrane potential is in-fluenced by the intracellular and extracellular relaxationtimes. Nonetheless, a more detailed study on how therelaxation times affect the virtual electrodes is needed tofully understand the hyperbolic bidomain model. Com-paring the stimulation patterns to experimental datacould even serve to calibrate the magnitude of the re-laxation times in the hyperbolic bidomain model. Despite significant progress in models of cardiac elec-trophysiology, the complex effects of nonlinearity andheterogeneity in the heart are far from being fully un-derstood. We have shown that the nonlinearities of theunderlying physics can give unexpected results, in con-tradiction to the linear case. Inductances in tissue prop-agation could play an important role in cardiac electricaldynamics, especially close to the wavefronts, where thefast currents responsible for the initiation of the actionpotential can give a small but non-negligible contribu-tion. The main phenomenon we observed in the hyper-bolic model is the enhancement of the conduction ve-locity at small relaxation times because of the presenceof the time derivative of the ionic currents. This phe-nomenon would allow electric signals to propagate at thesame speed with lower conductance as compared to thestandard models. Future experimental studies are neededto confirm the importance of these effects. ACKNOWLEDGEMENTS
The authors gratefully acknowledge research supportfrom NIH award HL117063 and from NSF award ACI1450327. The human atrial model used in this work waskindly provided by Prof. W. P. Segars. We also thankProf. C. S. Henriquez and Prof. J. P. Hummel for ex-tended discussions on atrial modeling. We also gratefullyacknowledge support by the libMesh developers in aid-ing in the development the code used in the simulationsreported herein.
Appendix A: Cable equation with inductances
We follow the derivation of the cable equation byKeener and Sneyd . As shown in Fig. 1, we assumethere exist inductance effects in both the extra- and in-tracellular axial currents. We have V e ( x ) − V e ( x + ∆ x ) = − L e ∂ t I e ∆ x − R e I e ∆ x, (A1) V i ( x ) − V i ( x + ∆ x ) = − L i ∂ t I i ∆ x − R i I i ∆ x, (A2)where I e and I i are the extracellular and intracellularaxial currents, respectively. Dividing by ∆ x and takingthe limit as ∆ x → , we find ∂ x V e = − L e ∂ t I e − R e I e , (A3) ∂ x V i = − L i ∂ t I i − R i I i . (A4)The minus sign on the right hand sides is a conventionthat ensures that positive charges flow from left to right.Applying Kirchhoff’s current law, we have that I e ( x ) − I e ( x + ∆ x ) + I t ∆ x = 0 , (A5) I i ( x ) − I e ( x + ∆ x ) − I t ∆ x = 0 , (A6)which, in the limit ∆ x → , becomes I t = ∂ x I e = − ∂ x I i . (A7)6Differentiating equations (A3) and (A4) with respect to x , and equation (A7) with respect to t , we find that forthe extracellular compartment, ∂ xx V e = − L e ∂ xt I e − R e ∂ x I e , (A8) ∂ xt I e = ∂ t I t , (A9)and for the intracellular compartment, ∂ xx V i = − L i ∂ xt I i − R i ∂ x I i , (A10) ∂ xt I e = − ∂ t I t . (A11)Substituting the mixed derivative of the extracellular andintracellular currents, we find ∂ xx V e = − L e ∂ t I t − R e I t , ∂ xx V i = L i ∂ t I t + R i I t . (A12)Defining V = V i − V e , and taking the difference betweenthe intracellular and the extracellular equations, we canwrite a single equation for V , ∂ xx V = ( L i + L e ) ∂ t I t + ( R i + R e ) I t . (A13)Recalling that I t = χ ( C m ∂ t V + I ion ) , (A14)we finally find the hyperbolic form of the cable equation, τ C m ∂ V∂t + C m ∂V∂t − ∂∂x (cid:18) D ∂V∂x (cid:19) = − I ion − τ ∂I ion ∂t , where we define the conductivity as D = 1 /χ ( R i + R e ) and the relaxation time by τ = L i + L e R i + R e . (A15) Appendix B: Exact solution for the piecewise-linear bistablemodel
Consider the simplified piecewise-linear bistablemodel, I ion ( V ) = kV − k [ V H ( V − V ) + V (1 − H ( V − V ))] . (B1)The nondimensional form of equation (B1) is the simpli-fied model proposed by McKean , ˆ I ion (cid:16) ˆ V (cid:17) = ˆ V − H (cid:16) ˆ V − α (cid:17) . (B2)Using (B1) in the hyperbolic monodomain model (21)and considering only one spatial dimension, we have µ ∂ ˆ V∂ ˆ t + (1 + µ ) ∂ ˆ V∂ ˆ t − ∂ ˆ V∂ ˆ x + ˆ V = 0 , ˆ V < α,µ ∂ ˆ V∂ ˆ t + (1 + µ ) ∂ ˆ V∂ ˆ t − ∂ ˆ V∂ ˆ x + ˆ V = 1 , ˆ V > α, (B3) where we have introduced the nondimensional variables ˆ t = 1 T ˆ t, ˆ x = 1 L x , ˆ V = V − V V − V . In particular, we define T = C m k , L = (cid:114) σkχ , µ = τ kC m , where µ is a nondimensional number that characterizesthe effect of the inductances in the system. Specifically, µ is the ratio between the relaxation time of the system andthe characteristic time of the reactions: the larger µ , themore important the effects of the inductances become.Introducing the change of coordinates z = x − ct, suchthat U ( z ) = U ( x − ct ) = ˆ V ( x, t ) , the system (B3) istransformed into (cid:40)(cid:0) c µ − (cid:1) U zz − c (1 + µ ) U z + U = 0 , z > , (cid:0) c µ − (cid:1) U zz − c (1 + µ ) U z + U = 1 , z < . (B4)Consider first the case z > . Using γ = c µ − and β = − c (1 + µ ) , this equation reads γV (cid:48)(cid:48) + βV (cid:48) + V = 0 , for which the roots of the characteristic polynomial are λ ± = − β ± (cid:112) β − γ γ , (B5)so that the solution is U ( z ) = A + e λ + z + A − e λ − z . It iseasy to verify that the case z < has solution U ( z ) = B + e λ + z + B − e λ − z + 1 . The global solution, therefore, is U ( z ) = (cid:40) A + e λ + z + A − e λ − z , z > ,B + e λ + z + B − e λ − z + 1 , z < . (B6)For (B6) to represent a propagating front, it is necessarythat the solution is real and bounded for any value of z .This implies that either λ − < < λ + or λ + < < λ − .Requiring β − γ > , we find the condition γ < β / ,which is always satisfied for the model (B4) because β − γ = c (1 − µ ) + 4 > . If c > then β < . Then λ − > and λ + < if c < (cid:114) µ . (B7)Imposing U ( ∞ ) = 0 , we find A − = 0 , and imposing U ( −∞ ) = 1 , we find B + = 0 , so that U ( z ) = (cid:40) A + e λ + z , z > ,B − e λ − z + 1 , z ≤ . (B8)To ensure the continuity of the solution at z = 0 , we fix U (0) = α , so that A + = α = 1 + B − and U ( z ) = (cid:40) αe λ + z , z > , ( α − e λ − z + 1 , z ≤ . (B9)7The first derivative of U is U (cid:48) ( z ) = (cid:40) αλ + e λ + z , z > , ( α − λ − e λ − z , z ≤ . (B10)Imposing the continuity of U (cid:48) at z = 0 , we find αλ + =( α − λ − and therefore γβ = (cid:0) α − α (cid:1) (2 α − . (B11)Using the definitions of γ and β in (B11), we finally ob-tain an expression for the speed c in terms of µ and α , c = (1 − α ) (cid:113) µ + ( α − α ) ( µ − ≤ (cid:114) µ . (B12)Notice that the characteristic propagation speed isbounded from above by (cid:112) /µ . Appendix C: Finite Element Discretization
The common practice to solve the coupled nonlinearsystem of cardiac electrophysiology is to split the mon-odomain (or bidomain) model from the ionic model. Mo-tivated by the work of Munteanu et al. , we will alsodecouple the two systems. For brevity, we show the dis-cretization only for the monodomain model with a first-order time integrator, because first-order time integratorsare still widely used in the community. The same ap-proach can be used to discretize the hyperbolic bidomainmodel. Given the subinterval [ t n , t n +1 ] , with t = 0 , wedefine two separate subproblems, one for the ionic modeland one for the monodomain equation connected by theinitial conditions. In the first step, we solve the ionicmodel (25) for w n +1 , and we compute the ionic currentsusing the updated values of the state variables. The ioniccurrents computed in this way are then used to solve themonodomain system (22)–(23).As shown in Spiteri and Dean , the most efficientnumerical algorithm to solve the stiff ODEs of the ionicmodel strongly depends on the model considered. We usethe forward Euler method for the simplified ionic modelsof Aliev-Panfilov and Fenton-Karma. For the ventricularionic model of ten Tusscher and Panfilov , we use theRush-Larsen method . The atrial ionic model of Grandiet al. has a stricter restriction on the time step size.Consequently, we use the Rush-Larsen method in com-bination to the backward Euler method for linear equa-tions. The time step size is defined as ∆ t = t n +1 − t n .To obtain a second-order time discretization scheme, theionic model can be solved using the explicit trapezoidalmethod (i.e., Heun’s method), which is a second-order ac-curate 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 finiteelement discretization. Because (22) is a simple ordinary differential equation, the computational costs for solving(21) or the system (22)–(23) is comparable. Denoting by ( v, w ) Ω = ´ Ω vw the L (Ω) inner product and using theboundary conditions (24), the Galerkin approximation ofequations (22)–(23) is to find V h , Q h ∈ S h such that (cid:18) ∂V h ∂t , ψ h (cid:19) Ω = (cid:0) Q h , ψ h (cid:1) Ω , (C1) (cid:18) τ C m ∂Q h ∂t + C m Q h , φ h (cid:19) Ω + (cid:0) D ∇ V, ∇ φ h (cid:1) Ω = (cid:0) I h ion , φ h (cid:1) Ω + (cid:18) τ ∂I h ion ∂t , φ h (cid:19) Ω , (C2)for all ψ h , φ h ∈ S h , where S h = (cid:8) v h ∈ C (cid:0) ¯Ω (cid:1) : v h (cid:12)(cid:12) K ∈ P ( K ) , ∀ K ∈ T h (cid:9) . (C3)In practice, we look for a continuous solution that is lin-ear in every simplex element K in the triangulation T h of Ω . Introducing the basis of nodal shape functions { N A } MA =1 , with M = dim ( S h ) , the solution fields arediscretized as V h = M (cid:88) A =1 N A V A , Q h = M (cid:88) A =1 N A Q A . Equations (C1) and (C2) yield the matrix system ˙ V = Q , (C4) C m M (cid:16) τ ˙ Q + Q (cid:17) + KV = F + τ L . (C5)We use the half-lumping scheme proposed by Path-manathan et al. , so that F = MI , L = MJ and ˙ V = Q , (C6) C m M L (cid:16) τ ˙ Q + Q (cid:17) + KV = M ( I + τ J ) , (C7)where [ M AB ] = ( N A , N B ) Ω , [ K AB ] = ( ∇ N A , D ∇ N B ) Ω , I and J are the nodal evaluation of I ion and ∂I ion /∂t , and M L is the lumped mass matrix obtained by row sum-mation of the mass matrix. Using a simple first-orderimplicit-explicit time integrator, the fully discrete systemof equations reads V n +1 = V n + ∆ t Q n +1 , (C8) C m ( τ + ∆ t ) M L Q n +1 + ∆ t KV n +1 = ∆ t M ( I ∗ + τ J ∗ )+ τ C m M L Q n , (C9)8where I ∗ and J ∗ indicate the values of the ionic currents,and of their derivatives, evaluated using the updated val-ues for w n +1 obtained by solving the ionic model. Intro-ducing (C8) in (C9), we arrive at a single system for Q n +1 , (cid:2) C m ( τ + ∆ t ) M L + ∆ t K (cid:3) Q n +1 = τ C m M L Q n − ∆ t KV n + ∆ t M ( I ∗ + τ J ∗ ) . (C10)The time derivative of the ionic currents, J ∗ , is approxi-mated nodally using the chain rule as ∂I ion ∂t ( t n ) ≈ ∂I ion ∂V (cid:0) V n , w n +1 (cid:1) Q n + ∂I ion ∂ w (cid:0) V n , w n +1 (cid:1) · g ( V n , w n +1 ) , (C11)where the function g ( V, w ) is the right hand side of theionic model system defined in equation (25). After solv-ing (C10), we update the voltage equation (C8).The scheme used to obtain (C8) and (C9) is the sim-plest method – ARS(1,1,1) – of a series of implicit-explicit(IMEX) Runge-Kutta (RK) algorithms that are popularfor hyperbolic systems . For a second-order time in-tegrator, we use the H-CN(2,2,2) scheme , which com-bines Heun’s method for the explicit part and the Crank-Nicholson method for the implicit part.We have not experienced any difference in time stepsize restriction between the parabolic and hyperbolicmodels. In fact, we have used an IMEX-RK methodwhere only the ionic currents and their time derivativesare treated explicitly in both parabolic and hyperbolicmodels. The time step size restriction is dictated only bythe ionic model. Nonetheless, the time step size shouldbe chosen to accurately capture the propagation of theaction potential. Denoting with v the conduction veloc-ity and with h m the smallest element size, we chose ourtime step size such that the CFL condition ∆ t ≤ h m /v holds.From a purely numerical perspective, in the first-orderscheme, the dissipative nature of the backward Eulermethod can be used to remove numerical oscillations,while the second-order Crank-Nicholson method will keepsuch oscillations. A dissipative second-order time inte-grator could be used in place of the Crank-Nicholsonmethod. In any case, the physical dissipation of the hy-perbolic monodomain model reduces such artifacts forfirst- and second-order schemes.One of the main difficulties arising from hyperbolic sys-tems is their tendencies to form discontinuities for whichnumerical approximations typically develop spurious os-cillations. Although the standard monodomain and bido-main models are parabolic systems, if the wavefront isnot accurately resolved then the upstroke of the actionpotential can act as a discontinuity. For this reason,the numerical methods used for cardiac electrophysiol-ogy should be already suitable for the hyperbolic systemsconsidered here. REFERENCES Abdusalam, H. A. 2004, ‘Analytic and approximate solutions forNagumo telegraph reaction diffusion equation’,
Applied Mathe-matics and Computation (2), 515–522. Ahrens, J., Geveci, B., Law, C., Hansen, C. and Johnson, C.2005, ‘36-paraview: An end-user tool for large-data visualiza-tion’,
The Visualization Handbook p. 717. Aliev, R. R. and Panfilov, A. V. 1996, ‘A simple two-variable model of cardiac excitation’,
Chaos, Solitons & Fractals (3), 293–301. Arthurs, C. J., Bishop, M. J. and Kay, D. 2012, ‘Efficient sim-ulation of cardiac electrical propagation using high order finiteelements’,
Journal of computational physics (10), 3946–3962. Ascher, U. M., Ruuth, S. J. and Spiteri, R. J. 1997, ‘Implicit-explicit runge-kutta methods for time-dependent partial differ-ential equations’,
Applied Numerical Mathematics (2-3), 151–167. Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune,P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D.,Kaushik, D., Knepley, M. G., McInnes, L. C., Rupp, K., Smith,B. F., Zampini, S. and Zhang, H. 2015 a , PETSc users manual,Technical Report ANL-95/11 - Revision 3.6, Argonne NationalLaboratory. URL: Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune,P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D.,Kaushik, D., Knepley, M. G., McInnes, L. C., Rupp, K., Smith,B. F., Zampini, S. and Zhang, H. 2015 b , ‘PETSc Web page’, . URL: Balay, S., Gropp, W. D., McInnes, L. C. and Smith, B. F. 1997,Efficient management of parallelism in object oriented numericalsoftware libraries, in E. Arge, A. M. Bruaset and H. P. Lang-tangen, eds, ‘Modern Software Tools in Scientific Computing’,Birkhäuser Press, pp. 163–202. Benson, A. P., Bernus, O., Dierckx, H., Gilbert, S. H., Green-wood, J. P., Holden, A. V., Mohee, K., Plein, S., Radjen-ovic, A., Ries, M. E. et al. 2010, ‘Construction and valida-tion of anisotropic and orthotropic ventricular geometries forquantitative predictive cardiac electrophysiology’,
Interface fo-cus p. rsfs20100005. Boscarino, S., Filbet, F. and Russo, G. 2016, ‘High order semi-implicit schemes for time dependent partial differential equa-tions’,
Journal of Scientific Computing (3), 975–1001. Boscarino, S. and Russo, G. 2009, ‘On a class of uniformly ac-curate imex runge–kutta schemes and applications to hyperbolicsystems with relaxation’,
SIAM Journal on Scientific Computing (3), 1926–1945. Bueno-Orovio, A., Cherry, E. M. and Fenton, F. H. 2008, ‘Min-imal model for human ventricular action potentials in tissue’,
Journal of Theoretical Biology (3), 544–560. Bueno-Orovio, A., Kay, D., Grau, V., Rodriguez, B. and Bur-rage, K. 2014, ‘Fractional diffusion models of cardiac elec-trical propagation: role of structural heterogeneity in disper-sion of repolarization’,
Journal of The Royal Society Interface (97), 20140352. Cattaneo, C. 1948, ‘Sulla conduzione del calore’,
Atti del Semi-nario matematico e fisico dell’Università di Modena , 83–101. Clapham, D. E. and Defelice, L. J. 1982, ‘Membrane BiolowSmall Signal Impedance of Heart Cell Membranes’,
The Jour-nal of Membrane Biology , 63–71. Clayton, R. H., Bernus, O., Cherry, E. M., Dierckx, H., Fenton,F. H., Mirabella, L., Panfilov, A. V., Sachse, F. B., Seemann, G.and Zhang, H. 2011, ‘Models of cardiac tissue electrophysiology:progress, challenges and open questions’,
Progress in biophysicsand molecular biology (1), 22–48. Colli Franzone, P., Pavarino, L. F. and Savaré, G. 2006, Compu-tational electrocardiology: mathematical and numerical model-ing, in ‘Complex Systems in Biomedicine’, Springer Milan, Mi- lano, pp. 187–241. Colli Franzone, P., Pavarino, L. F. and Scacchi, S. 2012, ‘Car-diac excitation mechanisms, wavefront dynamics and strength–interval curves predicted by 3d orthotropic bidomain simula-tions’,
Mathematical biosciences (1), 66–84. Colman, M. A. 2014, Atrial Fibrillation, pp. 201–225. DeHaan, R. L. and DeFelice, L. J. 1978, ‘Oscillatory Proper-ties and Excitability of the Heart Cell Membrane’,
TheoreticalChemistry: Periodicities in Chemistry and Biology , 181–233. Engelbrecht, J. 1981, ‘On theory of pulse transmission in a nervefiber’,
Proceeding Royal Sociaty London A , 195–209. Engelbrecht, J., Peets, T., Tamm, K., Laasmaa, M. and Vendelin,M. 2016, ‘On modelling of physical effects accompanying thepropagation of action potentials in nerve fibres’, arXiv preprintarXiv:1601.01867 . Fastl, T. E., Tobon-Gomez, C., Crozier, W. A., Whitaker, J.,Rajani, R., McCarthy, K. P., Sanchez-Quintana, D., Ho, S. Y.,O’Neill, M. D., Plank, G. et al. 2016, Personalized modelingpipeline for left atrial electromechanics, in ‘Computing in Cardi-ology Conference (CinC), 2016’, IEEE, pp. 225–228. Fenton, F. H., Cherry, E. M., Hastings, H. M. and Evans, S. J.2002, ‘Multiple mechanisms of spiral wave breakup in a model ofcardiac electrical activity.’,
Chaos (Woodbury, N.Y.) (3), 852–892. Fenton, F. H. and Karma, A. 1998, ‘Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Fila-ment instability and fibrillation’,
Chaos: An InterdisciplinaryJournal of Nonlinear Science (1), 20–47. Fort, J., Campos, D., González, J. R. and Velayos, J. 2004,‘Bounds for the propagation speed of combustion flames’,
Jour-nal of Physics A: Mathematical and General (29), 7185–7198. Fort, J. and Méndez, V. 2002 a , ‘Time-Delayed Spread of Virusesin Growing Plaques’, Physical Review Letters (17), 178101. Fort, J. and Méndez, V. 2002 b , ‘Wavefronts in time-delayedreaction-diffusion systems. Theory and comparison to experi-ment’, Reports on Progress in Physics . Gorecki, J. 1995, ‘Molecular dynamics simulations of a chemicalwave front’,
Physica D: Nonlinear Phenomena (1-2), 171–179. Grandi, E., Pandit, S. V., Voigt, N., Workman, A. J., Dobrev, D.,Jalife, J. and Bers, D. M. 2011, ‘Human Atrial Action Potentialand Ca2+ Model’,
Circulation Research . Griffith, B. E. and Peskin, C. S. 2013, ‘Electrophysiology’,
Com-munications on Pure and Applied Mathematics (12), 1837–1913. Hand, P. E. and Griffith, B. E. 2010, ‘Adaptive multiscale modelfor simulating cardiac conduction’,
Proceedings of the NationalAcademy of Sciences (33), 14603–14608. Hand, P. E., Griffith, B. E. and Peskin, C. S. 2009, ‘Derivingmacroscopic myocardial conductivities by homogenization of mi-croscopic models’,
Bulletin of mathematical biology (7), 1707–1726. Harrild, D. M. and Henriquez, C. S. 2000, ‘A computer modelof normal conduction in the human atria’,
Circulation research (7), e25–e36. Ho, S. Y., Anderson, R. H. and Sánchez-Quintana, D. 2002,‘Atrial structure and fibres: morphologic bases of atrial conduc-tion’,
Cardiovascular research (2), 325–336. Ho, S. Y., Cabrera, J. A. and Sanchez-Quintana, D. 2012, ‘Leftatrial anatomy revisited’,
Circulation: Arrhythmia and Electro-physiology (1), 220–228. Hodgkin, A. L. and Huxley, A. F. 1952, ‘A quantitative descrip-tion of membrane current and its application to conduction andexcitation in nerve’,
The Journal of physiology (4), 500. Hooks, D. A., Trew, M. L., Caldwell, B. J., Sands, G. B., LeGrice,I. J. and Smaill, B. H. 2007, ‘Laminar arrangement of ventricularmyocytes influences electrical behavior of the heart’,
CirculationResearch (10), e103–e112. Hurtado, D. E., Castro, S. and Gizzi, A. 2016, ‘Computationalmodeling of non-linear diffusion in cardiac electrophysiology: Anovel porous-medium approach’,
Computer Methods in Applied Mechanics and Engineering , 70–83. Jou, D., Casas-Vázquez, J. and Lebon, G. 1996, Extended Irre-versible Thermodynamics, in ‘Extended Irreversible Thermody-namics’, Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 41–74. Kaplan, S. and Trujillo, D. 1970, ‘Numerical Studies of the Par-tial Differential Equations Governing Nerve Impulse Conduction:The Effect of Lieberstein’s Inductance Term’,
Mathematical Bio-sciences , 379–404. Keener, J. P. and Sneyd, J. 2009,
Mathematical physiology ,Springer. King, J., Huang, C. and Fraser, J. 2013, ‘Determinants of my-ocardial conduction velocity: implications for arrhythmogenesis’,
Frontiers in Physiology , 154. Kirk, B. S., Peterson, J. W., Stogner, R. H. and Carey,G. F. 2006, ‘ libMesh : A C++ Library for Parallel AdaptiveMesh Refinement/Coarsening Simulations’,
Engineering withComputers (3–4), 237–254. http://dx.doi.org/10.1007/s00366-006-0049-3 . Knisley, S. B. 1995, ‘Transmembrane voltage changes duringunipolar stimulation of rabbit ventricle’,
Circulation Research (6), 1229–1239. Koch, C. 1984, ‘Cable theory in neurons with active, linearizedmembranes’,
Biological Cybernetics (1), 15–33. Krause, D., Dickopf, T., Potse, M. and Krause, R. 2015, ‘Towardsa large-scale scalable adaptive heart model using shallow treemeshes’,
Journal of Computational Physics , 79–94. Krishnamoorthi, S., Perotti, L. E., Borgstrom, N. P., Ajijola,O. A., Frid, A., Ponnaluri, A. V., Weiss, J. N., Qu, Z., Klug,W. S., Ennis, D. B. and Garfinkel, A. 2014, ‘Simulation Methodsand Validation Criteria for Modeling Cardiac Ventricular Elec-trophysiology’,
PLoS ONE (12), e114494. Krishnamoorthi, S., Sarkar, M. and Klug, W. S. 2013, ‘Numericalquadrature and operator splitting in finite element methods forcardiac electrophysiology’,
International Journal for NumericalMethods in Biomedical Engineering (11), 1243–1266. Krueger, M. W., Schmidt, V., Tobón, C., Weber, F. M., Lorenz,C., Keller, D. U., Barschdorf, H., Burdumy, M., Neher, P., Plank,G. et al. 2011, Modeling atrial fiber orientation in patient-specificgeometries: a semi-automatic rule-based approach, in ‘Interna-tional Conference on Functional Imaging and Modeling of theHeart’, Springer, pp. 223–232. Lemarchand, A. and Nawakowski, B. 1998, ‘Perturbation of localequilibrium by a chemical wave front’,
The Journal of chemicalphysics (16), 7028–7037. Lieberstein, H. M. 1967, ‘On the Hodgkin-Huxley partial differ-ential equation’,
Mathematical Biosciences (1), 45–69. Lieberstein, H. M. and Mahrous, M. A. 1970, ‘A source of largeinductance and concentrated moving magnetic fields on axons’,
Mathematical Biosciences (1-2), 41–60. MATLAB 2010, version 7.10.0 (R2010a) , The MathWorks Inc.,Natick, Massachusetts. McKean, H. P. 1970, ‘Nagumo’s equation’,
Advances in Mathe-matics (3), 209–223. Méndez, V. and Compte, A. 1998, ‘Wavefronts in bistable hyper-bolic reaction-diffusion systems’,
Physica A , 90–98. Méndez, V. and Llebot, J. E. 1997, ‘Hyperbolic reaction-diffusionequations for a forest fire model’,
Physical Review E (6), 6557. Munteanu, M., Pavarino, L. F. and Scacchi, S. 2009, ‘A scal-able newton–krylov–schwarz method for the bidomain reaction-diffusion system’,
SIAM Journal on Scientific Computing (5), 3861–3883. Nash, M. P. and Panfilov, A. V. 2004, ‘Electromechanicalmodel of excitable tissue to study reentrant cardiac arrhythmias’,
Progress in Biophysics and Molecular Biology (2), 501–522. Neu, J. C. and Krassowska, W. 1992, ‘Homogenization ofsyncytial tissues.’,
Critical reviews in biomedical engineering (2), 137–199. Niederer, S. A., Kerfoot, E., Benson, A. P., Bernabeu, M. O.,Bernus, O., Bradley, C., Cherry, E. M., Clayton, R., Fenton, F., Garny, A., Heidenreichm, E., Land, S., Maleckar, M., Path-manathan, P., Plank, G., Rodriguez, J. F., Roy, I., Sachse, F. B.,Seemann, G., Skavhaug, O. and Smith, N. P. 2011, ‘Verificationof cardiac tissue electrophysiology simulators using an N-versionbenchmark’,
Philosophical transactions. Series A, Mathematical,physical, and engineering sciences , 4331–4351. Papez, J. W. 1920, ‘Heart musculature of the atria’,
Develop-mental Dynamics (3), 255–285. Pashakhanloo, F., Herzka, D. A., Ashikaga, H., Mori, S., Gai,N., Bluemke, D. A., Trayanova, N. A. and McVeigh, E. R. 2016,‘Myofiber architecture of the human atria as revealed by submil-limeter diffusion tensor imaging’,
Circulation: Arrhythmia andElectrophysiology (4), e004133. Patelli, A. S., Dede, L., Lassila, T., Bartezzaghi, A. and Quar-teroni, A. 2017, ‘Isogeometric approximation of cardiac electro-physiology models on surfaces: An accuracy study with appli-cation to the human left atrium’,
Computer Methods in AppliedMechanics and Engineering , 248–273. Pathmanathan, P., Mirams, G. R., Southern, J. and Whiteley,J. P. 2011, ‘The significant effect of the choice of ionic currentintegration method in cardiac electro-physiological simulations’,
International Journal for Numerical Methods in Biomedical En-gineering (11), 1751–1770. Pezzuto, S., Hake, J. and Sundnes, J. 2016, ‘Space-discretizationerror analysis and stabilization schemes for conduction velocityin cardiac electrophysiology’,
International journal for numericalmethods in biomedical engineering . Potse, M., Dube, B., Richer, J., Vinet, A. and Gulrajani, R. M.2006, ‘A Comparison of Monodomain and Bidomain Reaction-Diffusion Models for Action Potential Propagation in the Hu-man Heart’,
IEEE Transactions on Biomedical Engineering (12), 2425–2435. Quarteroni, A., Lassila, T., Rossi, S. and Ruiz-Baier, R. 2017,‘Integrated Heart-Coupling multiscale and multiphysics modelsfor the simulation of the cardiac function’,
Computer Methods inApplied Mechanics and Engineering , 345–407. Rossi, S., Lassila, T., Ruiz-Baier, R., Sequeira, A. and Quar-teroni, A. 2014, ‘Thermodynamically consistent orthotropic ac-tivation model capturing ventricular systolic wall thickeningin cardiac electromechanics’,
European Journal of Mechanics -A/Solids , 129–142. Roth, B. J. 1997, ‘Electrical conductivity values used withthe bidomain model of cardiac tissue’,
IEEE Transactions onBiomedical Engineering (4), 326–328. Rush, S. and Larsen, H. 1978, ‘A practical algorithm for solvingdynamic membrane equations’,
IEEE Transactions on Biomedi-cal Engineering (4), 389–392. Sachse, F. B. 2004,
Computational Cardiology , Vol. 2966 of
Lec-ture Notes in Computer Science , Springer Berlin Heidelberg,Berlin, Heidelberg. Sánchez-Quintana, D., Pizarro, G., López-Mínguez, J. R., Ho,S. Y. and Cabrera, J. A. 2013, ‘Standardized review of atrialanatomy for cardiac electrophysiologists’,
Journal of cardiovas-cular translational research (2), 124–144. Scott, A. C. 1971, ‘Effect of the series inductance of a nerve axonupon its conduction velocity’,
Mathematical Biosciences (3-4), 277–290. Scott, A. C. 1972, ‘Transmission line equivalent for an unmyeli-nated nerve axon’,
Mathematical Biosciences (1-2), 47–54. Segars, W. P., Sturgeon, G., Mendonca, S., Grimes, J. and Tsui,B. M. W. 2010, ‘4d xcat phantom for multimodality imagingresearch’,
Medical physics (9), 4902–4915. Sepulveda, N. G., Roth, B. J. and Wikswo, J. P. 1989, ‘Currentinjection into a two-dimensional anisotropic bidomain’,
Biophys-ical journal (5), 987–999. Spach, M. S. 1983, ‘The discontinuous nature of electrical propa-gation in cardiac muscle - Consideration of a quantitative modelincorporating the membrane ionic properties and structural com-plexities the ALZA distinguished lecture’,
Annals of BiomedicalEngineering (3-4), 208–261. Spiteri, R. J. and Dean, R. C. 2010, ‘Stiffness Analysis of CardiacElectrophysiological Models’,
Annals of Biomedical Engineering (12), 3592–3604. Stan, D., del Teso, F. and Vázquez, J. L. 2014, ‘Finite and infinitespeed of propagation for porous medium equations with frac-tional pressure’,
Comptes Rendus Mathematique (2), 123–128. ten Tusscher, K. H. W. J. and Panfilov, A. V. 2006, ‘Alternansand spiral breakup in a human ventricular tissue model’, Amer-ican Journal of Physiology - Heart and Circulatory Physiology (3). Thomas, C. E. 1959, ‘The muscular architecture of the atria ofhog and dog hearts’,
Developmental Dynamics (2), 207–236. Tung, L. 1978, ‘A bi-domain model for describing ischemic my-ocardial d-c potentials.’. Vera, Z., Pride, H. P., Zipes, D. P. and Barr, R. C. 1995, ‘Reper-fusion arrhythmias: role of early afterdepolarizations studied bymonophasic action potential recordings in the intact canine heartduring autonomically denervated and stimulated states.’,
Jour-nal of cardiovascular electrophysiology (7), 532–43. Wikswo, J. P., Lin, S. F. and Abbas, R. A. 1995, ‘Virtual elec-trodes in cardiac tissue: a common mechanism for anodal andcathodal stimulation’,
Biophysical Journal (6), 2195–2210. Wikswo, J. P. and Roth, B. J. 2009, Virtual electrode theory ofpacing, in ‘Cardiac Bioelectric Therapy’, Springer, pp. 283–330. Zemskov, E. P., Tsyganov, M. A. and Horsthemke, W. 2015,‘Wavy fronts in a hyperbolic FitzHugh-Nagumo system and theeffects of cross diffusion’,
Physical Review E - Statistical, Non-linear, and Soft Matter Physics91