Yves Coudière
email: yves.coudiere@inria.fr

Myriam Rioux
email: myriam.rioux@mat.ulaval.ca

Virtual electrode polarization and current activation with monodomain equations

Keywords: virtual electrodes, monodomain model, anodal and cathodal stimulation, cardiac anisotropy

The bidomain model is nowadays one of the most accurate mathematical descriptions of the electrical activity in the heart. From now on, it was believed to be the only model for accurate simulations of cardiac muscle stimulation in a clinically relevant manner, i.e. through extracellular electrodes. In this paper, we develop a computationally efficient and accurate approximation of the bidomain model that allows for extracellular stimulation and accounts for unequal anisotropy ratios between intra-and extra-cellular media: the currentlifted monodomain model. We prove its use in the isolated heart by reproducing four types of extracellular activation that exhibit virtual electrode polarization. The simplicity of the code implementation and the fact that the computational cost is equal to the standard monodomain model henceforth give an excellent alternative to the bidomain model for expensive simulations like arrhythmia, fibrillation, and their respective treatments, e.g. signal guided catheter ablation and defibrillation.

Introduction

The electrical response of the heart tissue to a given external current stimulus remains a very active subject of research [START_REF] Trayanova | Modeling defibrillation of the heart: Approaches and insights[END_REF][START_REF] Tandri | Reversible cardiac conduction block and defibrillation with high-frequency electric field[END_REF][START_REF] Rantner | Three-dimensional mechanisms of increased vulnerability to electric shocks in myocardial infarction: Altered virtual electrode polarizations and conduction delay in the peri-infarct zone[END_REF][START_REF] Rantner | Terminating ventricular tachyarrhythmias using far-field low-voltage stimuli: Mechanisms and delivery protocols[END_REF]]. An electric discharge can as well induce ventricular fibrillation, or restore the normal sinus rhythm in a fibrillating heart [START_REF] Trayanova | Virtual electrodeinduced positive and negative graded responses: new insights into fibrillation induction and defibrillation[END_REF][START_REF] Trayanova | Bidomain model of defibrillation[END_REF]. Moreover, there is a plenty of observable phenomena when electric shocks are delivered, see [START_REF] Wikswo | Virtual electrode theory of pacing[END_REF] for a detailed review. For instance, experimental [START_REF] Efimov | Virtual electrode-induced phase singularity: a basic mechanism of defibrillation failure[END_REF][START_REF] Efimov | Direct evidence of the role of virtual electrode-induced phase singularity in success and failure of defibrillation[END_REF][START_REF] Sidorov | Examination of stimulation mechanism and strength-interval curve in cardiac tissue[END_REF][START_REF] Sidorov | Cathodal stimulation in the recovery phase of a propagating planar wave in the rabbit heart reveals four stimulation mechanisms[END_REF] and numerical studies [START_REF] Anderson | Termination of spiral waves with biphasic shocks[END_REF][START_REF] Franzone | Anode make and break excitation mechanisms and strength-interval curves: bidomain simulations in 3D rotational anisotropy[END_REF][START_REF] Franzone | Cardiac excitation mechanisms, wavefront dynamics and strength-interval curves predicted by 3D orthotropic bidomain simulations[END_REF] have demonstrated that extracellular cardiac stimulation by an unipolar electrode induces characteristic regions of membrane polarization other than those in the vicinity of the electrode. This whole class of phenomena is referred to as virtual electrode polarization (VEP). VEP have been hypothesized to be responsible for the failure of the shock to defibrillate [START_REF] Efimov | Virtual electrode-induced phase singularity: a basic mechanism of defibrillation failure[END_REF], a hypothesis that was supported with mathematical models [START_REF] Eason | Phase singularities and termination of spiral wave reentry[END_REF] and experimentally [START_REF] Efimov | Direct evidence of the role of virtual electrode-induced phase singularity in success and failure of defibrillation[END_REF].

The bidomain model is generally used for investigating defibrillation (e.g. [START_REF] Trayanova | Defibrillation of the heart: insights into mechanisms from modelling studies[END_REF]), where VEP may play an important role [START_REF] Trayanova | Virtual electrodeinduced positive and negative graded responses: new insights into fibrillation induction and defibrillation[END_REF]. Indeed, the bidomain representation of the cardiac tissue considers the superposition of two continuous media, the intracellular and the extracellular spaces. This allows electrical activation and shocks to be delivered in the extracellular space. It is also well known that myocardial conductivity properties are anisotropic, and furthermore, intraand extra-cellular media have different anisotropy ratios [START_REF] Guyader | Extracellular measurement of anisotropic bidomain myocardial conductivities. I. theoretical analysis[END_REF][START_REF] Roberts | Influence of cardiac fiber orientation on wavefront voltage, conduction velocity and tissue resistivity[END_REF][START_REF] Roberts | Effect of tissue anisotropy on extracellular potential fields in canine myocardium in situ[END_REF]. VEP is anisotropy related [START_REF] Sepulveda | Current injection into a twodimensional anisotropic bidomain[END_REF][START_REF] Wikswo | Virtual electrode theory of pacing[END_REF] and this condition is naturally accounted in the bidomain model.

The mathematical structure of the bidomain equations makes them computationally demanding. This is why the use of the bidomain approach is still debated for simulating real life situations of high complexity. For instance, modeling the effects of strong defibrillation shocks such as those used in electrical termination of fibrillations must take into account realistic fiber geometry, detailed membrane properties, microscopic tissue structure, inhomogeneity, and anisotropy of cardiac conductivities. These simulations describe phenomena that make sense only over many heart beats, and temporal resolution must be extremely fine to capture the ionic dynamics of the membrane. Moreover, meshes have to be extremely fine to capture the complex anatomical structure of the myocardium. For 2D and 3D simulations, the computations tend to be prohibitively expensive.

Much efforts are made to speed up the computations of the bidomain model. A first approach consists in using a range of more sophisticated one than the other numerical techniques. A good review of various strategies is given in [START_REF] Clayton | Models of cardiac tissue electrophysiology: Progress, challenges and open questions[END_REF]. A second approach is rather based on the choice of the models and their approximations. An example of approximation relies on the reduction of the bidomain model to one equation based on the concept of an activating function [START_REF] Sobie | A generalized activating function for predicting virtual electrodes in cardiac tissue[END_REF][START_REF] Skouibine | Increasing the computational efficiency of a bidomain model of defibrillation using a time-dependent activating function[END_REF]. Though it preserves most of bidomain properties and reduces the computation cost by a factor of two [START_REF] Skouibine | Increasing the computational efficiency of a bidomain model of defibrillation using a time-dependent activating function[END_REF], it remains unclear how the result could deviate from the bidomain solution [START_REF] Vigmond | Solvers for the cardiac bidomain equations[END_REF]. An other very common model that approaches the bidomain solution is the monodomain model, where no coupling with an extracellular domain is explicitly integrated. When the anisotropy ratios of the intra-and the extra-cellular media are considered equal and no current source is involved, both the monodomain and the bidomain solutions are the same. In general, monodomain and bidomain solutions are very similar [START_REF] Potse | A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF][START_REF] Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF]. An (not so significant) exception reported in the literature is about slight deviations in the trajectories of spiral wave tips [START_REF] Roth | Meandering of spiral waves in anisotropic cardiac tissue[END_REF], where no apparent difference in arrhythmia complexity is noticed. The monodomain is thus a well-accepted model for studying propagation phenomena/pathologies. In addition, it has been con-firmed by many authors that solving the monodomain model is about ten times faster than the bidomain model [START_REF] Bourgault | Comparing the bidomain and monodomain models in electro-cardiology through convergence analysis[END_REF][START_REF] Sundnes | On the computational complexity of the bidomain and the monodomain models of electrophysiology[END_REF][START_REF] Plank | Generation of histo-anatomically representative models of the individual heart: tools and application[END_REF]. For these reasons, monodomain model is widely used for investigating propagation phenomena on considerable periods of time, e.g. post-shock activity in arrhythmia and fibrillation studies [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation[END_REF][START_REF] Gharaviri | A computer model of endo-epicardial electrical dissociation and transmural conduction during atrial fibrillation[END_REF] on complex realistic domains [START_REF] Ten Tusscher | Organization of ventricular fibrillation in the human heart[END_REF][START_REF] Clayton | Vortex filament dynamics in computational models of ventricular fibrillation in the heart[END_REF][START_REF] Bishop | The role of fine-scale anatomical structure in the dynamics of reentry in computational models of the rabbit ventricles[END_REF]. We finally mention a hybrid approach inspired from both numerical techniques and choices of the models: a domain decomposition method [START_REF] Gerardo-Giorda | Optimized schwarz coupling of bidomain and monodomain models in electrocardiology[END_REF] where the monodomain equations are solved in subdomains where no current is applied, and the bidomain equations are solved elsewhere. In any of the approaches, there is always a compromise between the effort made to implement the numerical strategies and the validity of the model predictions.

Aside computational efficiency, studying the effect of a strong defibrillation shocks requires a model that accounts for biophysical realities. We summarize in three points these biophysical realities which have to be considered and modeled as is, and those which can be neglected.

1. Electrical activation and shocks are delivered in the extracellular space with leads that are usually placed on the endocardial heart surface. 2. The intra-and extra-cellular spaces are anisotropic with unequal anisotropy ratios, which directly causes the VEP. 3. The effect of unequal anisotropy ratio can be neglected in absence of external electrical source.

The extracellular stimulation (1.) as well as the unequal anisotropy ratio (2.) are naturally accounted in the bidomain formulation. Meanwhile, the standard monodomain model accounts for none of these. As a consequence, the bidomain model is generally considered inescapable for investigating defibrillation, while the monodomain model remains an attractive alternative because of its low computational cost.

In this paper, we propose to derive a monodomain model in the presence of intra-or extra-cellular current stimulations. The main idea relies on the linearity of the quasi-static balance between the intra-and extra-cellular spaces in the bidomain equations. Following the superposition principle, we split its solution into two parts: the excitation part, directly related to the external sources, and the propagation parts, related to the propagation dynamics in the tissue. The excitation part is kept unchanged and thus accounts properly for extracellular stimulations and unequal anisotropy ratios (see (1.) and (2.) above), while the propagation part is replaced by a monodomain approximation (following (3.) above).

We end up with a monodomain equation with a new source term that is obtained by solving the electrical balance for the external stimulation only. This amounts to lift the external current stimulation before solving the monodomain equation, hence the model is called current-lifted (CL) monodomain model. Naturally, this additional source term has to be computed only when the stimulation is active. Moreover, it can be precomputed because it depends only on the tissue, the lead and stimulation protocol. This preprocessing is negligible compared to the time step iterations. The CL-monodomain model is finally a computationally efficient alternative to the bidomain model for long-duration and clinically relevant numerical simulations involving any kind of electrical stimulation involving VEP.

The details of the mathematical derivation are given in section 2. In section 3, we present the computational setup and the numerical methods. We finally prove the use of the CL-monodomain model in section 4, where the four well known VEP phenomena, specifically anode/cathode make/break are observed when applying extracellular current with an endocardial surface electrode.

Models equations for an isolated tissue with stimulation electrodes

We assume that the excitable cardiac tissue occupies an open, bounded, and connected subset Ω of R 3 . For sake of simplicity, we assume that it is insulated and we denote by Γ its boundary (figure 1).

Following the classical bidomain modeling approach, the heart system is described by the intra-and extra-cellular potentials u i (t, x) and u e (t, x), the electrophysiological ionic state variables collected in a vector w(t, x), all defined in Ω. As usual, the transmembrane potential is defined by the difference v = u i -u e in Ω.

We assume that some external electrodes are used for the stimulation of the cardiac tissue. These electrodes are located usually on the boundary Γ, and acts as surface sources of current with known densities related either to the intra-or the extra-cellular media, respectively denoted by j i (t, x) and j e (t, x), defined on Γ. For completeness, we also consider volume sources of currents with densities s i (t, x) and s e (t, x), respectively associated to the intra-and extracellular media. They are defined in Ω.

Derivation of the current-lifted monodomain equations

Here we will observe that the stimulation by electrodes can be exactly written as an integral source term in the bidomain equation. By linearity, the main part of the bidomain equations can be replaced by its monodomain approximation while retaining the integral source term, leading to a monodomain equation that accounts for the virtual electrodes phenomena through the integral source term.

The bidomain equations

The bidomain model for the heart Ω with all these sources consists in the following system of PDEs, for x ∈ Ω and t > 0,

χ (C m ∂ t v + I ion (v, w)) -∇ • (σ i ∇u i ) = s i , (1)
χ (C m ∂ t v + I ion (v, w)) + ∇ • (σ e ∇u e ) = -s e , (2)
∂ t w + g(v, w) = 0, (3)
completed by the boundary conditions

σ i ∇u i • n = j i , σ e ∇u e • n = j e , (4)
on Γ and for t > 0. For compatibility reasons, we require that, for all t > 0,

Ω (s i + s e ) dx + Γ (j i + j e ) dγ(x) = 0. (5)
The initial condition is given by

v(0, x) = v 0 (x), w(0, x) = w 0 (x) in Ω. (6)
The data are the current sources j i , j e , s i , and s e . The parameters of the equations are the ratio of surface of cell membrane per unit volume of tissue χ, the cell membrane capacitance per unit area C m , and the anisotropic intra-and extra-cellular electrical conductivity matrices σ i (x) and σ e (x). The electrophysiological behavior of the system is characterized by the functions I ion (v, w) and g(v, w) given for instance by the Beeler-Reuter [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibers[END_REF] or the Mitchell-Schaeffer [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF] models. The myocardium is assumed to have a laminar structure described by the data of a local orthonormal basis (a 1 (x), a 2 (x), a 3 (x)) where a 1 (x) is along the fiber and (a 2 (x), a 3 (x)) is along the laminae. Hence the electrical conductivities are explicitly given by

σ i,e (x) = 3 d=1 σ d i,e a d (x) T a d (x),
where σ 1 i,e ≥ σ 2 i,e ≥ σ 3 i,e > 0 are given scalar electrical conductivities (see table 1).

The equations (1) to (3) are called the parabolic-parabolic formulation of the bidomain model. Here, we will use the equivalent parabolic-elliptic formulation below (obtained with u i = v + u e ), for x ∈ Ω and t > 0:

χ (C m ∂ t v + I ion (v, w)) -∇ • (σ i ∇u i ) = s i , (7)
-∇ • (σ i ∇u i ) -∇ • (σ e ∇u e ) = s i + s e , (8)
∂ t w + g(v, w) = 0, (9)
From the variational form of the bidomain equations based on this parabolicelliptic formulation, and using the boundary conditions ( 4)), we get the variational formulation, for a.e. t > 0:

χ C m d dt (v, φ) + (I ion (v, w), φ) + a i (v + u e , φ) = (s i , φ) + Γ j i φ (10)
a i (v + u e , φ e ) + a e (u e , φ e ) = (s i + s e , φ e ) + Γ (j i + j e )φ e (11)
d dt (w, ψ) + (g(v, w), ψ) = 0, (12)
for all test functions φ, φ e ∈ H 1 (Ω), and ψ ∈ L 2 (Ω). Here (•, •) denotes de L 2 (Ω) inner product and a i , a e are bilinear forms defined on H 1 (Ω) by a i,e (u, φ) = Ω σ i,e ∇u • ∇φ. The functional spaces are chosen according to [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF].

The current-lifted monodomain equations

We apply the superposition principle in order to split equation [START_REF] Sidorov | Cathodal stimulation in the recovery phase of a propagating planar wave in the rabbit heart reveals four stimulation mechanisms[END_REF] in two equations. We search for the solution u e into the form u e = u e solutions of the following two equations: a i (u (1) e , φ e ) + a e (u (1) e , φ e ) = -a i (v, φ e ),

a i (u (2) e , φ e ) + a e (u (2) e , φ e ) = (s i + s e , φ e ) + Γ

(j i + j e )φ e . (14)
Remark that the potential u

(1) e

depends only on the transmembrane potential v, and is computed exactly like for the bidomain system of equations without external stimulation, while the potential u

(2)
e depends only on the current sources but not on the transmembrane potential v.

Consequently, we call the potential u

"transmembrane-dependent extracellular potential", and the potential u e in equation ( 10) leads to the following equation for the transmembrane potential v:

χ C m d dt (v, φ) + (I ion (v, w), φ) + a i (v + u (1) e , φ) = (s i , φ) -a i (u (2) e , φ) + Γ j i φ. (15)
Together with the equations ( 13) and ( 14) on u

(1)

e , u (2)
e , and equation ( 9) on w, it is an equivalent statement of the complete bidomain equations ( 7) to [START_REF] Efimov | Direct evidence of the role of virtual electrode-induced phase singularity in success and failure of defibrillation[END_REF].

In order to obtain the "current-lifted (CL) monodomain equations", we replace the bidomain operator defined by Av, φ := a i (v+u

(1) e , φ) for u (1)
e solution of ( 13) by its monodomain approximation, defined by

A m v, φ = a m (v, φ) := Ω σ m ∇u • ∇φ with σ m = (σ i ) -1 + (σ e ) -1 -1 .
For equal anisotropy ratio, σ i = λσ e for some fixed λ > 0, these operators coincide and σ m = 1 λ+1 σ i = λ λ+1 σ e . In general, both operators propagate planar waves at the same velocity in the eigendirections a i (i = 1, 2, 3). The choice of σ m can be adjusted to fit other propagation properties, see [START_REF] Coudière | Optimal monodomain approximations of the bidomain equations used in cardiac electrophysiology[END_REF].

The CL-monodomain equations read:

χ C m d dt (v, φ) + (I ion (v, w), φ) + a m (v, φ) = (s i , φ) -a i (u (2) e , φ) + Γ j i φ (16)
d dt (w, ψ) + (g(v, w), ψ) = 0, (17)
a i (u (2) e , φ e ) + a e (u (2) e , φ e ) = (s i + s e , φ e ) + Γ

(j i + j e )φ e , (18)
where the current-lifted extracellular potential u

(2) e is the same as the one calculated in the original bidomain equations ( 13), [START_REF] Franzone | Cardiac excitation mechanisms, wavefront dynamics and strength-interval curves predicted by 3D orthotropic bidomain simulations[END_REF], and [START_REF] Eason | Phase singularities and termination of spiral wave reentry[END_REF].

Note that equations ( 16) and ( 17) are exactly the weak formulation of the usual monodomain equations, with an additional source term depending only on u

(2) e . These equations entirely decouple from the equation ( 14) that defines of u

(2) e .

The strong form of the current-lifted monodomain equations

It can be easily seen that the corresponding strong formulation of this model reads:

χ (C m ∂ t v + I ion (v, w)) = ∇ • (σ m ∇v) + s i + ∇ • σ i ∇u (2) e , (19)
∂ t w + g(v, w) = 0, (20)
for x ∈ Ω and t > 0, with the boundary condition (for x ∈ Γ and t > 0)

σ m ∇v • n = j i -σ i ∇u (2) e • n. (21)
In these equations, the current-lifted extracellular potential is computed by solving independently the elliptic problem below, for all t > 0:

-∇ • (σ i + σ e ) ∇u (2) e = s i + s e , in Ω, (22)
(σ i + σ e ) ∇u (2) e • n = j i + j e , on Γ. ( 23
)
The solution u

(2) e depends only on the applied current s i , s e , j i , and j e , in a linear way. Remark that the current-lifted extracellular potential u [START_REF] Tandri | Reversible cardiac conduction block and defibrillation with high-frequency electric field[END_REF] e is defined up to an additional constant, for each time t > 0 . Anyway, the final source term ∇ • σ i ∇u

(2) e in equation [START_REF] Roberts | Effect of tissue anisotropy on extracellular potential fields in canine myocardium in situ[END_REF] does not depend on this constant. If necessary, the constant can be fixed for instance by grounding u

(2) e at some fixed location x 0 ∈ Ω, so that u

(2) e (t, x 0 ) = 0 for all t > 0. If no external stimulation is applied (s i = s e = 0 and j i = j i = 0), then the problem has the trivial solution u

(2) e = 0, and the current-lifted model ( 19)- [START_REF] Sepulveda | Current injection into a twodimensional anisotropic bidomain[END_REF] degenerates to the usual monodomain model.

The computation of u

(2) e is required only during the external stimulation, which has a negligible effect on the overall computation.

Furthermore, for practical reasons, only surface currents are usually applied, and only through the extracellular media, which specifically means that s i = s e = 0 and j i = 0. If the electrodes are fixed, it is not restrictive to assume that j e (t, x) = H(t) je (x), where je (x) is a normalized surface density of current, and H(t) is the temporal modulation of the current, i.e. defined by the stimulation protocol. For instance, in the setup of section 3 (figure 1), we consider two extracellular electrodes placed on the endo and epicardial surfaces Γ endo and Γ ref . A total current I tot is applied through these electrodes that uniformly distribute the current on their surface. In this case, the normalized density of current is

je (x) =      Itot |Γ endo | for x ∈ Γ endo , -Itot |Γ ref | for x ∈ Γ ref , 0 elsewhere on Γ. (24)
Instead of ( 22)-( 23) defined for all t > 0, the time independant normalized problem

-∇ • (σ i + σ e ) ∇ū (2) e = 0, in Ω, (25)
(σ i + σ e ) ∇ū (2) e • n = je , on Γ, (26)
are solved once and the current-lifted extracellular potential is finally computed explicitly as

u (2) e (t, x) = H(t)ū (2) e (x). (27)

Reconstruction of the total extracellular potential

For practical reason, it can be important to reconstruct the complete extracellular potential from the CL-monodomain model. It is given by u e = u e can be recovered by solving the electrostatic balance equation [START_REF] Franzone | Anode make and break excitation mechanisms and strength-interval curves: bidomain simulations in 3D rotational anisotropy[END_REF]. The corresponding strong form of this problem reads, for all time t > 0:

-∇ • (σ i + σ e )∇u (1) e = ∇ • (σ i ∇v) in Ω, (28)
(σ i + σ e )∇u (1) e • n = -σ i ∇v • n on Γ. ( 29
)
As the potential u

e , the solution u

(1) e (t, x) is defined up to an additional constant at each time t > 0. Similarly, it can be fixed by grounding the potential u [START_REF] Trayanova | Modeling defibrillation of the heart: Approaches and insights[END_REF] e at the given location x 0 ∈ Ω, and we also have u

-∇ • (1 + λ)σ e ∇u (1) e = ∇ • (λσ e ∇v)
in Ω,

(1 + λ)σ e ∇u (1) e • n = -λσ e ∇v • n on Γ, (30)
for all time t > 0, or equivalently -∇ • (σ e ∇Ψ(t)) = 0 in Ω with the boundary condition σ e ∇Ψ(t) • n = 0 on Γ, where

Ψ(t) = (1 + λ)u (1)
e (t) + λv(t). This shows that Ψ(t) is a constant, independent of x ∈ Ω.

Again, this constant in space Ψ(t) can be fixed by grounding the extracellular potential u

(1) e at the fixed location x 0 ∈ Ω, meaning that ∀x ∈ Ω, (1 + λ)u (1) e (t, x) + λv(t, x) = (1 + λ)u (1) e (t, x 0 ) + λv(t, x 0 ) = λv(t, x 0 ).

Finally, the transmembrane extracellular potential is

∀t > 0, ∀x ∈ Ω, u (1) e (t, x) = - λ 1 + λ (v(t, x) -v(t, x 0 )) . (32)

Resolution steps

The complete resolution of the CL-monodomain can be carried out following the three decoupled steps summarized below. The lifting operation makes independent any calculation implying currents (step 1), that is performed once before all time iterations of step 2. The CL-monodomain model is then not significantly more computationally expensive than the standard monodomain model.

More complex situations are allowed easily. For instance consider several electrodes with normalized densities je,k (x) (k = 1, 2 . . .), and different time activation sequences H k (t). An elemental contribution ū(2) e,k is then computed for each electrode, and the total current-lifted extracellular potential is the linear combination of theses contributions u

e (t, x) = k H k (t)ū (2)
e,k (x).

Step 1 Preprocessing the CL-monodomain model Require: The data of the normalized source of current je (x) and timecourse H(t). ū(2) e (x) ← the solution of equation ( 25) and [START_REF] Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF]. return u

(2)

e (t, x) ← H(t)ū (2)
e (x), the current-lifted extracellular potential.

Step 2 Processing the CL-monodomain model

Require:

The current-lifted extracellular potential u Iterate in time to solve the monodomain equation ( 19) and [START_REF] Sepulveda | Current injection into a twodimensional anisotropic bidomain[END_REF] with boundary condition [START_REF] Clayton | Models of cardiac tissue electrophysiology: Progress, challenges and open questions[END_REF].

return The CL-monodomain solution v(t, x) and w(t, x).

Step

Postprocessing the CL-monodomain model [optional]

Require:

The current-lifted extracellular potential u

e (t, x). Require: The transmembrane voltage v(t, x). Require: The location x 0 of the ground for extracellular measures.

if Equal anisotropy ratio assumption, σ i = λσ e then u

(1)

e ← -λ 1+λ (v(t, x) -v(t, x 0 )), else u (1)
e ← the solution of equation ( 28) and [START_REF] Sundnes | On the computational complexity of the bidomain and the monodomain models of electrophysiology[END_REF].

end if return u e ← u (1)
e + u

e , the total extracellular potential. We consider also a 2D geometry, which is a cross section (z = 0) of the 3D domain. This represents a transversal cross section of the heart ventricular walls. These geometries are inspired from the modeling framework described in [START_REF] Franzone | Cardiac excitation mechanisms, wavefront dynamics and strength-interval curves predicted by 3D orthotropic bidomain simulations[END_REF].

Experimental setup

Model parameters

We use the Mitchell-Schaeffer (MS) model [START_REF] Tolkacheva | Analysis of the Fenton-Karma model through an approximation by a one-dimensional map[END_REF], a two-variable model describing the ionic activity of a cell with the transmembrane potential v and a scalar recovery variable w. Referring to equations ( 19) and ( 20), the MS source terms writes [START_REF] Boulakia | Towards the numerical simulation of electrocardiograms[END_REF][START_REF] Relan | Quantitative comparison of two cardiac electrophysiology models using personalisation to optical and MR data[END_REF][START_REF] Relan | Parameter estimation of a 3D cardiac electrophysiology model including the restitution curve using optical and MR data[END_REF][START_REF] Marchesseau | Personalization of a cardiac electromechanical model using reduced order unscented kalman filtering from regional volumes[END_REF]. It has several qualities that make it a candidate of choice for matching several action potential features [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF] (conduction speed, time scales, restitution). However, we understand it is not suitable for replicating realistic responses to defibrillation shocks and we refer the reader to a good review [START_REF] Trayanova | Defibrillation of the heart: insights into mechanisms from modelling studies[END_REF] giving arguments based on experimental observations. The MS model is finally a good compromise for computation speed while allowing us to achieve our goal: to prove the relevance and accuracy of the CL-monodomain model.

I ion (v, w) = 1 τ in wv 2 (v -1) + 1 τ out v, (33)
g(v, w) = 1 τopen (w -1) for v < v gate , 1 τ close w for v ≥ v gate . ( 34
)
The conductivities eigenvalues σ 1,2,3

i,e are reported in table 1 and all the other model parameters are reported in table 2. Remark that the first set of conductivity eigenvalues is that of [START_REF] Franzone | Cardiac excitation mechanisms, wavefront dynamics and strength-interval curves predicted by 3D orthotropic bidomain simulations[END_REF] but with σ 2

i,e = σ 3 i,e . In the 2D geometry, the longitudinal and normal fiber directions are along x and y respectively. In 3D, direction x is along the fibers and directions y and z are across fibers.

Virtual electrodes stimulation protocol

In this section, we consider surface electrodes only. This choice is made by simplicity but other kind of sources may be considered as well.

For make stimulation the tissue is initially at rest (v = 0, w = 0.99). A super-threshold total current I tot is turned on at t = 0 for a duration of 10 ms. This means that the function H in equation ( 24) is H(t) = 1 for 0 ≤ t ≤ 10 and 0 afterwards. An excitation wavefront starts propagating during the excitation, and the propagation continues at the turnoff of the electrode.

For break stimulation the tissue is initially refractory (v = 0, w = 0.1), and the current I tot is turned on until a relative refractory state is estab-lished. Throughout the application of the current, the recovery variable w slowly changes. In the vicinity of the endocardial electrode develop excited regions that are kept refractory through time. Anywhere else becomes potentially excitable because of the slow increase of the recovery variable w. No propagation of an excitation wavefront is possible at that step, as the excited regions are refractory. A sort of steady state is then observed for the transmembrane potential v. At the turnoff of the electrode the excited regions return to the resting state as they are refractory. The excitable surrounding tissue then depolarizes because of diffusion, and an excitation wavefront sweeps the cardiac domain.

Numerical methods

Discretization method

We consider unstructured finite element meshes of Ω which characteristics are specified in table 3 for the 2D and 3D computations. The bidomain and monodomain weak formulations are discretized on these meshes by the P1-Lagrange finite element technique, and the second-order completely implicit Gear time integration scheme. The first time step is calculated with a Crank-Nicholson time integration scheme. We choose a fixed time-step of ∆t = 0.1 ms, coherent with the dynamics of the MS model and the spatial resolution of the diffusion. The resulting linear systems are solved with Newton's method where 3-4 iterations are performed each time step for a stopping criteria of 1 • 10 -6 . The result is a sequence of discrete solutions of the form (v n h ) n≥0 that approximate the transmembrane potential v(t n , •) at times t n = n∆t. All results were obtained with the finite element library MEF++ [START_REF]Projet MEF++[END_REF].

Isochrons of depolarization

Depolarization times t depol (x) are recorded in order to compare the results of the simulations. For each position x ∈ Ω, they provide the time when the transmembrane potential v(t, x) crosses a threshold v th = 0.13 from below, i.e. from the rest state to the excited state. Given x ∈ Ω, for computation of depolarization times, we use the Lagrange polynomial of degree 2 that interpolates the three points (t n-1 , v n-1 h (x)), (t n , v n h (x)), and (t n+1 , v n+1 h (x)):

p 2 (t) = v h [t n-1 ] + v h [t n-1 , t n ](t -t n-1 ) + v h [t n-1 , t n , t n+1 ](t -t n-1 )(t -t n ).
Here v h [. . .] denotes the Newtons divided-differences formula coefficients associate to the interpolation points. Remark that Gear time integration scheme, which is used in our simulations, is derived with the same polynomial. At each time step, we look for a root of p 2 (t) -v th in the interval [t n-1 , t n ] and for which p ′ 2 (t) > 0 (v is increasing). The value of t depol assigned on every degree of freedom is the first root encountered during the iterations.

Results

In this section, we validate the use of the CL-monodomain model by reproducing experimental phenomena that are initiated by extracellular cathodal or anodal stimulation of cardiac tissue. We first review in section 4.1.1 the solution of the standard monodomain and bidomain models with the activation by a transmembrane stimulus. With the use of the bidomain model, we review in section 4.1.2 the effect of the equal/unequal anisotropy ratios between the intraand the extra-cellular media.

The CL-monodomain model is used thereafter to prove the existence of the virtual cathode and anode when applying extracellular current with an endocardial electrode. The four well known VEP phenomena, i.e. anode/cathode make/break, are reproduced in section 4.2 with the CL-monodomain solution and compared with the bidomain solution.

Comparison of models under current activation

The usual way to activate the monodomain model is to apply a transmembrane current. We first evaluate the effect of the application of a transmembrane current on a 2D Cartesian slab of myocardium and show that no VEP phenomenon is possible even with the bidomain model. We then review the effect of equal/unequal anisotropy ratio when an extracellular current is applied with the bidomain model on the same 2D geometry.

Activation by a transmembrane stimulus

The monodomain equation with the MS ionic model and a transmembrane current applied on the surface writes, for x ∈ Ω and t > 0,

χ (C m ∂ t v + I ion (v, w)) -∇ • (σ m ∇v) = 0, (35)
∂ t w + g(v, w) = 0, (36)
with the boundary condition written on v,

σ m ∇v • n = j s , (37)
for x ∈ Γ and t > 0. Here j s is the surface density of the applied current. A transmembrane current does not result from a contribution of the intra or extracellular space. It is rather a source of current at the membrane, which is fictitious. In the bidomain model such a current must leave the intracellular space as it enters the extracellular space, and vice versa. In this case, the boundary condition (4) for the bidomain becomes:

σ e ∇u e • n = j s , σ i ∇u i • n = -j s .
which is conformal with the compatibility condition.

The fiber structure laminae are generally tangent to the boundary where the electrodes are located (Γ endo and Γ ref ). The third eigendirection a 3 of σ 3 i,e is then aligned with the normal n to Γ. It is also the third eigendirection of σ m = σ i -1 + σ e -1 -1 . In this case, we can write σ i n = σ 3 i n, σ e n = σ 3 e n, and

σ m n = σ 3 m n, where σ 3 m = σ 3 i -1 + σ 3 e -1 -1 = σ 3 i σ 3 e σ 3 i +σ 3 e
. Hence, the bidomain boundary conditions (4) are σ 3 i ∂ n u i = j s and σ 3 e ∂ n u e = -j s and a straightforward computation shows that

σ 3 i ∂ n u i = j s σ 3 e ∂ n u e = -j s ⇒ σ 3 i σ 3 e ∂ n u i = σ 3 e j s σ 3 i σ 3 e ∂ n u e = -σ 3 i j s ⇒ σ 3 i σ 3 e σ 3 i + σ 3 e ∂ n (u i -u e ) = j s . ( 38
)
This is exactly the monodomain boundary condition [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibers[END_REF].

Hence we compare the monodomain model with the applied density of current j s to the bidomain with the same applied density of current j s . We use a current with a total magnitude I tot = 0.00625 mA during 10 ms. This current is applied through the two electrodes Γ endo and Γ ref from figure 1 so that the density of current j s is also given by equation [START_REF] Vigmond | Solvers for the cardiac bidomain equations[END_REF]. The epicaridal and endocardial electrodes then apply a current of -0.00625 mA/cm and 0.125 mA/cm respectively.

The solutions are compared in figure 2. The isochrons of depolarization are almost superposed. The depolarization waves are almost planar and propagate in the second eigendirection of the fiber structure. The error of propagation in this case should be nearly 0, as shown in [START_REF] Coudière | Optimal monodomain approximations of the bidomain equations used in cardiac electrophysiology[END_REF]. As a matter of fact, the observed relative error on depolarization times is at most 0.025%. It is also coherent with previous observations of differences between the monodomain and bidomain models [START_REF] Potse | A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF][START_REF] Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF].

Effect of unequal anisotropy ratios

In this section, the current I tot = 0.03125 mA is delivered for a duration of 10 ms. According to equation ( 24), for such a positive current I tot > 0, the anodal electrode is Γ endo . For the bidomain model with the equal anisotropy ratio, the eigenvalues conductivities are those specified in the second row of table 2.

Figure 3 illustrates clearly the important role of considering unequal anisotropy ratios between the intra-and the extra-cellular spaces. With unequal anisotropy ratios, two depolarized regions form during the stimulation on both sides of the anode (A) in the longitudinal direction. These regions are called the virtual cathodes (VC). The depolarization starts up from the virtual cathodes and also near the reference cathode (C) at the epicardial surface. The solution of the CLmonodomain model is very close to the solution of the bidomain model. Indeed, we observe that the formation of virtual cathodes and the propagation patterns are very similar. For the equal anisotropy setup, no VEP effect is observed, and the activation starts uniformly from the reference cathode.

VEP phenomena with the CL-monodomain model

The following results compare the solutions of the CL-monodomain solution with the bidomain solutions on the 3D geometry. The initial conditions and the stimulation protocol are summarized in table 4 for all the four VEP, where only the total current is given. The endo and epicardial density of currents are given in equation [START_REF] Vigmond | Solvers for the cardiac bidomain equations[END_REF]. The same stimulation protocol is applied for both the bidomain and the CL-monodomain simulations and the results are compared. VEP Initial conditions (for x ∈ Ω) Simulation protocol CM v 0 = 0, w 0 = 0.99 I tot (t ≤ 10 ms) = -0.125 mA AM v 0 = 0, w 0 = 0.99 I tot (t ≤ 10 ms) = 0.125 mA CB v 0 = 0, w 0 = 0.1 I tot (t ≤ 100 ms) = -2.5 mA AB v 0 = 0, w 0 = 0.1 I tot (t ≤ 100 ms) = 2.5 mA 4.2.1. Cathode Make (CM) and Anode Make (AM) A current stimulus is applied at the myocardial surface and is then turned off. For the CM simulations, the transmembrane potential u at 10 ms is shown in figures 4 and 5 for, respectively, the bidomain and the CL-monodomain models. Both the bidomain and the CL-monodomain solution show the typical pattern of excitation: a dog-bone shaped virtual cathode (depolarized region) in the tangential direction and two adjacent virtual anodes (hyperpolarized regions) in the longitudinal direction. Activation times are shown in figure 6. For the AM simulations, the transmembrane potential u of the bidomain model is shown in figure 7 at 10 ms, i.e. just before the turnoff of the electrode. The transmembrane potential u of the CL-monodomain model is very similar and will not be shown to avoid redundancy. Again, both the bidomain and the CL-monodomain solution show the typical pattern of excitation: two virtual cathodes in the fiber longitudinal direction. Activation times are shown in figure 8. Very similar propagation patterns are observed for both models.

Cathode Break (CB) and Anode Break (AB)

For break stimulations, the tissue is initially refractory. A strong stimulus is applied and a relative refractory state rapidly establishes. the excitable zone for the bidomain and the CL-monodomain simulation respectively. The refractory zone is defined by {x ∈ Ω; w(t, x) < 0.2267} according to the remark from section 3.1.2. It covers completely the excited region; hence, no AP propagation is possible. The excited region returns to the resting state at the turnoff of the electrodes as it is refractory. The excitable surrounding tissue then depolarizes because of diffusion, and an excitation wavefront sweeps the cardiac domain.

Activation thresholds

Lifting or not the current does not change the solution, as the error between the CL-monodomain and the bidomain models comes from the approximation of the propagation term only. However, the activation threshold may be affected by propagation and boundary effects. In this section, we attempt to find the activation thresholds of an extracellular stimulation of 10 ms in the 2D geometry by varying the magnitude of the current I tot . Activation thresholds for the AM and CM are presented in table 5. When the stimulation is made by the epicardial electrode (e.g. first row in table 5), it initiates a nearly plane wave propagating in an eigendirection (y) of the fiber arrangement. It can be proved that in homogeneous (but anisotropic) media, plane wave monodomain solutions propagating in the conductivity eigendirections are also bidomain solutions, as long as the harmonic average of the intra-and extra-cellular conductivities is used in the monodomain conductivity [START_REF] Coudière | Optimal monodomain approximations of the bidomain equations used in cardiac electrophysiology[END_REF]. Hence for cathodal epicardial stimulation (first row), the CL-monodomain solution is so close to the bidomain solution that we measure the same activation threshold. However, this is not the case for stimulation through the virtual cathodes or the endocardial cathode, where the stimulating wavefronts are curved. The difference between the activation threshold of the bidomain and the CL-monodomain models is actually greater when comparing virtual cathode stimulation (of about 5%) or cathodal endocardial stimulation (of about 9%). This could be explained by diffusion ar-guments. Indeed, for finding the activation threshold we are dealing with very low currents, but just high enough so that the diffusion does not extinguish the activation. Estimating the diffusion term with a monodomain approximation may impact the magnitude of the diffusion and so the activation threshold.

Discussion and conclusion

Discussion

Our model consists in the monodomain equation with the additional source term ∇ • σ i ∇u , where the so-called current-lifted extracellular potential u (2) e is precomputed. This term depends only on the electrode of stimulation itself, and on the structural properties of the tissues (specifically, the unequal anisotropy ratio), but not on its electrophysiology. We can determine the response of the tissue to an extracellular current source located anywhere in the tissue, then by superposition we can determine the total response by summing these responses weighted by the time varying magnitude of the current. Without exception u

(2) e originates from no approximation at all. Hence this component presents all aspects of the extracellular bidomain stimulation, e.g the VEP effects related with unequal anisotropy ratios in the intra-and extra-cellular spaces.

We verified first the preliminary statements of the derivation:

• lifting or not the current in the bidomain model gives the same solution up to the machine precision, and

• unequal anisotropy ratios are necessary to the apparition of VEP.

We then confronted the CL-monodomain model with the four well-known extracellular activation procedures: cathode make, anode make, cathode break, and anode break. The activation patterns associated with the VEP effects were found to be very similar as those of the bidomain model. The activation times were also in good agreement, the variation being attributable to the different propagation of the activating wavefront. We finally denoted slight discrepancies in bidomain and CL-monodomain thresholds activation, again related to the propagation.

Limitations. This was not in the scope of this paper to investigate the propagation discrepancies between the bidomain and the monodomain solutions. This issue has been thoroughly studied for instance in [START_REF] Potse | A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF][START_REF] Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF][START_REF] Bourgault | Comparing the bidomain and monodomain models in electro-cardiology through convergence analysis[END_REF]. In general very minor differences are measured. The range of error is usually considered as negligible for such a field of application (biological sciences) where the tolerance on the data errors is much larger [START_REF] Bourgault | Comparing the bidomain and monodomain models in electro-cardiology through convergence analysis[END_REF]. This is why the monodomain model continue to be preferred for high resolution simulations over extended periods of time [START_REF] Ten Tusscher | Organization of ventricular fibrillation in the human heart[END_REF][START_REF] Clayton | Vortex filament dynamics in computational models of ventricular fibrillation in the heart[END_REF][START_REF] Bishop | The role of fine-scale anatomical structure in the dynamics of reentry in computational models of the rabbit ventricles[END_REF][START_REF] Gharaviri | A computer model of endo-epicardial electrical dissociation and transmural conduction during atrial fibrillation[END_REF]. In this study, we have disregarded the presence of external passive regions representing extracardiac media like torso, lungs and blood in the heart cavities.

However, the derivation for the fully coupled intra/extra-cardiac model is also possible. The main difficulty relies on transmission conditions that can affect the wavefront curvature and trajectory near the boundary. These phenomena are referred to as bath-loading effects. A heuristic augmented monodomain approach has been developed in [START_REF] Bishop | Representing cardiac bidomain bath-loading effects by an augmented monodomain approach: application to complex ventricular models[END_REF] to address this issue quite efficiently.

Conclusion

The goal of this study was to derive in a rigorous theoretical framework a new model that provides a computationally efficient alternative to the bidomain model in the isolated heart. The resulting CL-monodomain model allows clinically relevant stimulation through the extracellular space, and accurately accounts for unequal anisotropy ratio related effects like VEP.

In consequence, the CL-monodomain model extends the reach of modeling applications of the monodomain approach to the investigation of complex mechanisms of stimulation and termination, e.g. VEP due to fiber curvature, VEP effects that may induce reentrant waves [START_REF] Sambelashvili | Dynamics of virtual electrode-induced scroll-wave reentry in a 3D bidomain model[END_REF] and are known to be a possible cause of arrhythmia and fibrillation, success or failure of defibrillation [START_REF] Anderson | Termination of spiral waves with biphasic shocks[END_REF], post-shock activation patterns, graded responses [START_REF] Trayanova | Defibrillation of the heart: insights into mechanisms from modelling studies[END_REF], etc.

( 2 )

 2 e "current-lifted extracellular potential". Substituting u e by u

e

  is already known. In general the transmembrane-dependent extracellular potential u[START_REF] Trayanova | Modeling defibrillation of the heart: Approaches and insights[END_REF]

( 1 )

 1 e (t, x 0 ) = 0 for all t > 0. The bidomain and monodomain equations without source of current are equivalent under the equal anisotropy ratio assumption. Recovering u (1) e can thus be performed using the same assumption. If σ i = λσ e with λ the anisotropy ratio, then the u (1) e -recovery problem simplifies into

3. 1 .

 1 Modeling setup 3.1.1. Computational domainsWe consider first a 3D domain of dimensions 1×1×0.45 cm 3 for representing a Cartesian slab of the heart ventricular walls (see figure1). All the stimulation pulses are delivered with extracellular surface currents j e through a Neumann boundary condition. A fixed current of I tot (in mA) is applied between two electrodes as depicted on figure1. The smaller electrode of dimension 0.05×0.05 cm 2 injects current on Γ endo , i.e. at the center of the endocardial interface (upper boundary). The larger electrode represents a reference electrode and covers the whole epicardial surface Γ ref (lower boundary). This modeling setup is made to approach an invasive activation procedure for stimulating the ventricle of a patient in diagnosis or treatment, e.g. during a signal guided catheter ablation, where the catheter electrode activates the endocardium by injecting an extracellular current.

Figure 1 :

 1 Figure 1: 3D (left) and 2D (right) computational domains. Surface electrodes are in red.

Figure 2 :

 2 Figure 2: Colormap: Depolarization times (ms) for the bidomain model with the application of a total transmembrane current of 0.0625 mA during 10 ms. Contours: isochrons of depolarization for the bidomain and the CL-monodomain models, respectively in black and white, and almost superposed.

Figure 3 :

 3 Figure 3: Colormap: depolarization times (ms) of the bidomain model with the application of a total extracellular current of 0.03125 mA during 10 ms. Contours: isochrons of depolarization for the bidomain model (black) and the CL-monodomain model (white). Left: Unequal anisotropy ratios. Right: Equal anisotropy ratio.

Figure 4 :

 4 Figure 4: Cathode make with the bidomain model: Transmembrane potential u at 10 ms, i.e. just before the turnoff of the electrode, for three cross sections of the 3D domain (left, top right and bottom right): the endocardial surface, the longitudinal and the transversal fiber directions.

Figure 5 :

 5 Figure 5: Cathode make with the CL-monodomain model: Transmembrane potential u at 10 ms for the same three cross sections as in figure 4.

Figures 9 (

 9 CB) and 12 (AB) show the transmembrane potential u at 100 ms, i.e. just before the turnoff of the stimulation. The black and white lines separate the refractory zone and

Figure 6 :

 6 Figure 6: Cathode make: Depolarization times (ms) for the bidomain model (colormap), and isochrons for the bidomain (black contours) and CL-monodomain (white contours) models.

Figure 7 :

 7 Figure 7: Anode make with the bidomain model: Transmembrane potential u at 10 ms.

Figure 8 :

 8 Figure 8: Anode make: Depolarization times (ms) for the bidomain model (colormap), and isochrons for the bidomain (black contours) and CL-monodomain (white contours) models.

  Figures 10 (CB) and 13 (AB) shows how the activation starts near the electrode from the edges of the refractory zones and propagates through the excitable regions. Activation times are shown in figures 11 and 14. Both the bidomain and the CL-monodomain solutions show the typical pattern of break excitation.

Figure 9 :

 9 Figure 9: Cathode break: Transmembrane potential u of the bidomain model after a cathodal activation of -2.5 mA during 100 ms. The black and white lines separate the refractory zone (R) and the excitable zone (E) for, respectively, the bidomain and the CL-monodomain models.

Figure 10 :

 10 Figure 10: Cathode break: Evolution of the transmembrane potential u for the CLmonodomain model after the turnoff of the cathode. The solution on the endocardial surface is shown (from left to right, top to bottom) 2, 4, 6, 8, 10, 25 ms after the turnoff.

Figure 11 :

 11 Figure 11: Cathode break: Depolarization times (ms) for the bidomain model with isochrons (bidomain: black contours, CL-monodomain: white contours). The densely packed contours indicates that the depolarization wavefront propagates around the refractory zone.

Figure 12 :

 12 Figure 12: Anode break: Transmembrane potential u of the bidomain model after an anodal activation of 2.5 mA during 100 ms. The black and white lines separate the refractory zone (R) and the excitable zone (E) for, respectively, the bidomain and the CL-monodomain models.

Figure 13 :

 13 Figure 13: Anode break: Evolution of the transmembrane potential u for the CL-monodomain model after the turnoff of the anode. The solution on the endocardial surface is shown (from left to right, top to bottom) 2, 4, 6, 8, 10, 25 ms after the turnoff.

Figure 14 :

 14 Figure 14: Anode break: Depolarization times (ms) for the bidomain model with isochrons (bidomain: black contours, CL-monodomain: white contours). The densely packed contours indicates that the depolarization wavefront propagates around the refractory zone.

Table 1 :

 1 Conductivities [S/m].

		σ 1 i	σ 2 i	σ 3 i	σ 1 e	σ 2 e	σ 3 e
	Unequal anisotropy ratio 0.23 0.0244 0.024 0.15 0.10	0.10
	Anisotropy ratio	-	9.5	9.5	-	1.4	1.4
	Equal anisotropy ratio	0.23 0.024 0.024 0.15 0.016 0.016
	Anisotropy ratio	-	9.5	9.5	-	9.5	9.5

Table 2 :

 2 Model parameters.

	C m	χ	τ in	τ out	τ open τ close v gate
	[F/m 2 ] [m -1 ]	[ms]	[ms]	[ms]	[ms]	[-]

1 • 10 -2 2 • 10 5 0.315 5.556 94.94 168.5 0.13

In the MS model, the rest state is at (v, w) = (0, 1) and the excited state is around v = 1, or very close to the curve v = 1/2 + 1/2 1 -4τ in /(τ out w). The refractory state is determined by the value of the recovery variable w, specifically when it is below 4τ in /τ out . Values of table 2 indicates that the refractory state is given by w < 0.2267.

Many authors have used the MS model for realistic clinical applications

Table 3 :

 3 Experimental meshes

	# Elements # Vertices Maximal Mesh size
			mesh size near the electrode
	3D 334, 400	60, 315	0.25 mm	0.025 mm
	2D 20, 985	10, 732	0.01 mm	0.001 mm

Table 4 :

 4 VEP stimulation protocol.

Table 5 :

 5 Activation thresholds.

	Endocardial Type of	Bidomain CL-monodomain
	electrode	activation site [mA]	[mA]
	Anode	C	0.00459	0.00459
		VC	0.00509	0.00488
	Cathode	C	0.000883	0.000800

Acknowledgements

This work was partially supported by an ANR grant part of "Investissements d'Avenir" program reference ANR-10-IAHU-04 and by an ANR grant within the program "Modèles Numériques 2013" reference ANR-13-MONU-0004-01. We also thank the financial support from FRQNT and from the GIREF for the computational resources and inestimable help from research professionals.