Yvan Kraepiel

Marie-Anne Barny
email: barny@upmc.fr

Gram-negative phytopathogenic bacteria, all hemibiotrophs after all?

Traditionally, the lifestyles of plant pathogens have been divided into distinct groups predicated on nutrient acquisition and the viability of host tissue. Biotrophs extract nutrients from living cells while necrotrophs feed off killed cells.Necrotrophy is defined as growth and nutrition of the pathogen on dead or dying plant material.

Hemibiotrophsshare characteristics with both biotrophs and necrotrophs, initially invading cells that are maintained alive prior to a transition to a necrotrophic life style in which nutrients are obtained from killing host cells.

The biotrophic/necrotrophic terminology is not adapted to bacteria Althoughthis terminology has been developed with pathogenic fungi, it is sometimes used to describe bacterial life styles. However, when applied to bacteria, the picture is not so clear, and we believe that the biotrophic/necrotrophic terminology should not be used for bacteria. A brief survey of the literature shows that Pseudomonas syringae is often refereed to as biotrophic or hemibiotrophic but is also occasionally described as partly necrotrophic or even necrotrophic. Xanthomonas spp. are often referred to as biotrophic. Erwiniaamylovora is often referred to as necrogenic, which cautiously describes its capacity to induce necrotic symptoms regardless of its life style, but it could also be described as necrotrophic or hemibiotrophic. The soft rot pathogens, Pectobacterium spp. and Dickeya spp. are mostly described as necrotrophs. Moreover, both biotrophy and necrotrophy have been assigned to Ralstoniasolanacearum in the literature. The delivered message is obviously rather confusing with bacteria.Why is this?

To understand, let us look back to the initial descriptions of biotrophic and necrotrophic life style. These contrasted lifestyles were initially described following microscopic observation of plant colonisation by pathogenic fungi (Lo Prestiet al, 2015). As paradigm ofbiotrophic fungus, the maize pathogen Ustilagomaydisinitiates colonisation by an initial intracellular growth, during which its hyphae is encased by the plant plasma membrane, and later switches to predominantly intercellular growth. In addition, the biotrophic tomato pathogen Cladosporiumfulvum remains exclusively in the extracellular compartment of tomato leaves. During the colonisation processes by these two biotrophic fungi, the plant cells remain alive. These microscopic observations were the firstto indicate that plant pathogens may possess tools to avoid plant cell death.Howeverhyphaes of necrotrophic fungi such as Botrytis cinerea or Sclerotiniasclerotiorum, after a short initial period of intercellular subcuticulargrowth, kill epidermal cells as soon as they penetrate inside the plant cells.In this mode of infection, the death of host plant cells precedes or accompanies colonization by the pathogen. Hemibiotrophic fungi such as Colletotrichumspp. or Magnaportheoryzae initially develop an intracellular bulged hypha encased by the plant plasma membranewhich does not kill the plant cell, and later switch to a thin intracellular necrotrophic hypha. When applied to bacteria, these assignments were not driven by close microscopic examination of the disease processes. Indeed, gram negative plant pathogenic bacteria remain extracellular throughout the infection process and it is not possible to classified them as biotrophs or necrotrophs on the basis of microscopic examination, as it has been performed with fungi. This could explain the difficulties to clearly assign a biotrophic or necrotrophic life style to bacteria. One bacterial exception that is not discussed in this opinion letter is Agrobacterium which could be considered as a true biotroph because it does not kill its host plant cells to proliferate, but rather induces the development of plant tumours producing metabolites that it can catabolize.

Why has this confusing terminology been applied to bacteria?

The explanationprobably liesin the fact thatit is important for plant pathologists to be able to distinguish between biotrophic and necrotrophic pathogens,asplant defences against biotrophs and necrotrophsare distinct and antagonistic. Because biotrophic pathogens require a living host, localized controlled cell deathtriggered through specific recognition of the invading pathogen by specific plant resistance (R) proteins forms part of an effective defence strategy. Necrotrophic pathogens, however,actively kill host tissue and therefore programmed cell death initiated by the plant is intuitively not an effective strategy to limit necrotrophic pathogen growth. It is also widely accepted that the salicylicacid (SA)-dependent defence signalling is effective against biotrophic pathogens, whereasjasmonicacid(JA) andethylene(ET)-dependent defence responses are efficient against necrotrophic pathogens. Determining the life style of a given plant pathogen is therefore important. This explains why plant pathologistshave attempted to classify bacteria as necrotrophs and biotrophs, as it has been performed with fungi.

However, when applied to bacteria, the biotrophic/necrotrophic terminology is not always consistent with experimental results, leading scientists to struggle in endless discussions to fit their data with the dogma. For example,itmaybestated in an article dealing with the influence on plant defence of the ERF96 transcription factor belonging to the ethylene responsive factor family« 'overexpression of ERF96 increased Arabidopsis resistance to necrotrophic pathogens such as the fungi Botritiscinerea and Pectobacteriumcarotovorum ssp. carotovorumbacteria.

However, Arabidopsis overexpressing ERF96 was more sensitive to hemibiotrophic bacteria P.syringaepv.

Tomato DC3000 » [START_REF] Catinot | ETHYLENE RESPONSE FACTOR 96 positivelyregulatesArabidopsisresistance to necrotrophicpathogens by direct binding to GCC elements of jasmonate -and ethylene-responsive defencegenes[END_REF] necrotrophAlternariabrassicicola ». In our opinion, the efficiency of a given plant defense mechanism against a given pathogen should not dictate the terminology used to describe the pathogen's life style, this only add confusion. Indeed, JA-related defence mechanisms, which include the production of most major classes of secondary metabolites and defence-related proteins, trigger efficient defences against many pathogens, whatever their life style [START_REF] Campos | Jasmonate-Triggered Plant Immunity[END_REF]. Furthermore, the dichotomy JA-necrotroph/SA-biotroph has mainly emerged from studies with Arabidopsis whereas studies with other species such as Monocots illustrate a much more contrasted reality [START_REF] De | Hormone defense networking in rice: tales from a different world[END_REF]. Even on Arabidopsis, the picture is not simple for bacterial pathogens. For example, JA-related defences are partly efficient during the symptomatic macerating phase induced by the soft rot pectinolitic bacteria D. dadantii, as expected for a "necrotrophic" pathogen. However plant necrosis, which is supposed to benefit to necrotrophic pathogens is also able to block efficiently D. dadantiiat the onset of infection (Kraepielet al., 2011). This clearly indicates that D. dadantii can't be classified so easily as biotroph or necrotroph. Similarly both SA and JA/ET pathways are efficient against Pectobacteria and the relative importance of both types of defence depends on the stage of the infection (Davidssonet al. 2013).

Biotrophic and hemibiotrophic type III-dependentphytobacteria, in addition to defence suppression, alsoinducecelldeath

Thenecrotrophic/ biotrophicclassification terminology is often superimposed on the main weapons deployed bydifferent bacterial pathogen.One one hand bacterial pathogensthat rely on a functional type III secretion system (T3SS) deploy a battery of injectedtype III effectors(T3Es)mostly involved in the suppression of plant defences for the benefit of the pathogen. This underscores the importance of the initial biotrophic development of these bacteria and explains why these bacteria are mostly described as biotrophic or hemibiotrophic. On the other hand, soft rot pectinolytic bacteria, which secrete a large set of plant cell wall-degrading enzyme (PCWDEs) through the type II secretion system (T2SS), are particularly effective in macerating the host tissues and in obtaining nutrients from the dead cells,explaining why they have often been referred to as necrotrophic bacteria. However, pathogenicity is not that simple and, as described below,a closer look at the pathogenicweaponsdeployed by bacterialpathogensrules out this oversimplified terminology.

Let us examine first the pathogenicity of Pseudomonas, Xanthomonas, Ralstonia and Erwiniaspecies thatrely on a functional T3SS. This T3SS allows the injection of T3Es inside the plant host cells.Asinjected T3Es can also be detected by plant R proteins, pathogens, during evolution, either loose the recognized effector or acquire new effectors to avoid R proteinsrecognition.This has beenconceptualized in the wellknownzig-zag model of plant pathogen evolution. As a consequence, biotrophs pathogens often deploy a large set of T3Es, are recognized by many plant R proteinsand have a narrow host range. Already,we can see that the pathogenicity of E. amylovora and R. solanacearumdoes not fit totally to this scheme. E.amylovora has a narrow host range, a characteristic of biotrophs, but there are no known resistance genes against this bacterium, a characteristic mostly shared by necrotrophs [START_REF] Malnoy | Fireblight: appliedgenomic insights of the pathogen and host[END_REF]. R.

solanacearum injects a plethora of T3Es involved in the suppression of plant defences, a characteristic of biotrophs, and yet this bacterium has a very large host range similar necrotrophs (Pouyemiro and Genin, 2009).Although the zig-zag model is extremely powerful in deciphering the mechanisms by which plant-pathogenic bacteria suppress plant defences, this paradigm has also limited our ability to see beyond it. Indeed, biotrophic and hemibiotrophic type III-dependent phytobacteria, in addition to defence suppression, also induce cell death. The need for an efficient mode of nutrient acquisition is a possible trigger for the transition to necrotrophy, and cell death is generally the ultimate response to pathogen attack. Type III-dependent phytobacteria harbour in their genomes several weapons to achieve this goal.First,some T3Es are involved in cell death induction during disease development. Interestingly, often, these cell death eliciting T3Es are key players for pathogen growth in planta, suggesting that induction of cell death is important for nutrients acquisition during infection. For example, PthA, which belongs to the transcription activator-like effectors (TALE) T3E family, is essential for pathogenicity of X .citri. PthA is involved in programmed killing of host cells and ectopic expression of pthA in citrus cells is sufficient to cause typical disease symptoms: division, enlargement and death of host cells. PthA is required for the production of necrotic cankers on all species of citrus attacked by X. citri. TALEs are DNA-binding proteins with a modular DNA-binding domain. The DNA binding domain is predictable, which simplifies elucidation of TALE function in planta, and it was shown that PthA activates the expression of the transcription factor CsLOB1, which coordinates pustule formation [START_REF] Boch | TAL effectors-pathogenstrategies and plant resistance engineering[END_REF]. Another T3Es family involved in cell death elicitation is the AvrE-like T3E family. This family is widespread among type IIIdependent phytobacteria and plays a crucial role for bacterial growth in planta.

Interestingly, T3Es of this family play a dual role during the disease process; they inhibit SA-mediated plant defences and interfere with vesicular trafficking, but also elicit electrolyte leakage and a slow plant cell death that participatesin nutriment release [START_REF] Degrave | The AvrE superfamily: ancestral type III effectors involvedin suppression of pathogen-associated molecular pattern-triggered immunity[END_REF]. These dual effects suggest that AvrE-like effectors may be involved in the transition from biotrophy to necrotrophy.

In addition to T3Es, type III-dependent phytobacteria harbour other virulence factors devoted to plant cell death. Pseudomonas, Xanthomonas, and Ralstonia species harbour T2SS and PCWDEs. These PCWDEs have often been suggested to participate mainly in the saprophytic life of these pathogens. However, although mutations affecting the T2SS of P. syringae have not been reported, R. solanacearum, X. campestrispv. campestris, X.

oryzaepv. oryzae and X. campestrispvvesicatoriaT2SS mutants are affected in virulence, indicating that T2SS is also an important virulence factor for these bacteria [START_REF] Szczesny | Functionalcharacterization of the Xcs and Xps type II secretionsystemsfrom the plant pathogenicbacteriumXanthomonascampestrispvvesicatoriaNew[END_REF].

Other factors involved in plant cell death are the toxins produced by P. syringaestrains.

The coronatine toxin has been shown to act as a mimic leading to the activation of JAdependent defences and therefore to the suppression of antagonistic SA-dependent defences. Moreover, this toxin is known to injure plant cells and to be involved in disease development and symptoms. For example, the coronatine-defective P.

syringaepv. atropurpurea mutant strain does not cause any symptoms,whereas the coronatine-producing strain induces a severe water-soaking symptoms, on oat [START_REF] Yao | ApoptoticCellDeathis a Common Response to PathogenAttack in Oats[END_REF]. In addition, others pathovars of P. syringae produce other toxins, such as tabtoxin, phaseolotoxin or mangotoxin, which can interfere with the nitrogen metabolism of the host, causing amino acid deficiencies in host cells and concomitant accumulation of nitrogen-containing intermediates that can be metabolized by the pathogen as a nitrogen source. This metabolic imbalance of amino acids leads to chlorosis and even necrosis symptoms in the host plant, and probably aids pathogen growth, because of the release of nutrients [START_REF] Arrebola | Chemical and Metabolic Aspects of AntimetaboliteToxinsProduced by Pseudomonas syringaePathovarsToxins[END_REF].

In pectobacterium and Dickeya soft rot necrotrophicpathogens, in addition to brute force killing, the asymptomaticbiotrophic phase plays an important role Now Let us examine the situation of pathogenic soft rot enterobacteria of genera Pectobacteriumand Dickeya.For these bacteria, as described above, pathogenicity relies mostly on a T2SS whichallows secretion of a large arsenal of PCWDEs,such as pectinases, cellulases, hemicellulases, and proteases. Once expressed, these large amounts of PCWDEs rapidly disrupt host cell integrity and promote rotting. However soft rot enterobacteriacould also be found in latent asymptomatic infections on many host crops. Theselatent infections, which arehighly dependent on environmental conditions, can last for severalmonths with the pathogen reaching high population levels in the absence of visible disease symptoms [START_REF] Liu | Quorum SensingCoordinates Brute Force and Stealth Modes of Infection in the Plant PathogenPectobacteriumatrosepticum[END_REF]. These asymptomatic infections play an important role in the disease dynamicsof natural infections. Because these asymptomatic infections are difficult to reproduce in laboratory conditions, the exact nature of thisbiotrophic phase remains elusive.One of the key elements to explain the switch between this asymptomatic biotrophic phase and the soft-rotting necrotrophic phase is probably the fine-tuning of PCWDEs production.Indeed, the production of PCWDEs is strictly controlled in a population density-dependent manner through quorum sensing (QS) regulationin Pectobacterium, whereas the PecS global regulatorpreventsthe premature expression of D. dadantiiPCWDEsat the beginningof theinfectionprocessbeforethe appearance of symptoms [START_REF] Mhedbi-Hajri | PecS is an important player in the regulatory network governing the coordinated expression of virulence genes during the interaction between Dickeyadadantii 3937 and plants[END_REF]. This fine regulation is probably necessary to prevent premature activation of plant defences, as the action of PCWDEs releases cell wall fragments which trigger defence responses in the host plant [START_REF] Davidsson | Pathogenicity of and plant immunity to soft rot pectobacteria[END_REF].

In adition to PCWDEs regulation, quorum sensing in Pectobacteriumalso regulates othervirulence determinants that could be involved in plant defence suppression.Among them, the T3SS and the AvrE-like T3EDspE/Ahave been found to be important for P.

atrosepticumvirulence. As effectors of this family are involved in the suppression of SAmediated plant defence, it would be interesting to test whether DspE/A proteins plays a similar role in P. atrosepticum. Similarly, the P. atrosepticumtype VI secretion system (T6SS)is also regulated through quorum sensing and T6SS mutants are also affected in virulence. Moreover, P. wasabiaeharborstwo T6SS machineries thathave been shownexperimentally to have partiallyoverlappingfunctionsduringpotato infection [START_REF] Davidsson | Pathogenicity of and plant immunity to soft rot pectobacteria[END_REF]. The T6SS, like T3SS, allows the injection of type VI effectors proteins (T6Es) into eukaryotic cells. The exact nature and role of the T6Es injected by P.

atrosepticum or P. wasabiaeare still unknown, but it can be speculated that they may be involved in the suppression of plant defences.Other traits identified through genomic studies could also benefit the bacterium. Among them, a type IV secretion system (T4SS)

and a putative polyketidephytotoxin (encoded by the cfacluster) have been shown to contribute to the virulence of P.atrosepticum. Once again, thenatureand the role ofthematerialtranslocatedthroughtheT4SSofP.atrosepticuminside the plant cellremainto be determined.

The diverse life style of phytopathogenic bacteria should be seen as a continuum of hemibiotrophic pathogens.

Irrespective of the bacteria considered, after an initial biotrophic phase, in which the bacteria avoid the elicitation of or suppress plant defences, a switch to necrotrophic stage,in which the bacteria actively kill plant cells, is observed. Both phases are required to induce disease although the length of time of each phase will vary between bacteria.

The length of each phase will also vary for the same bacteria between experimental conditions, as the inoculum concentration, age of the plant and the inoculation procedure will influence the outcome of the interaction. Another difficulty arises because the switch between the initial biotrophic phase and the later necrotrophicphase cannot be visualised through microscopic examination, as it has been performed with fungi. This makes the necrotroph/biotroph classification highly hazardous for phytopathogenic bacteria. Furthermore, all phytopathogenic bacteria possess weapons to supress or avoid premature elicitation of plants defences, indicating they initially develop on living hosts, and weapons involved in cell death induction,probably resulting in nutrient acquisition. Therefore, phytopathogenic bacteria should always be considered as hemibiotrophs. Although the biotrophic phase has been thoroughly described for type III-dependent phytobacteria, it has been overlooked for type IIdependent soft rot pectinolytic bacteria, probably because the asymptomaticbiotrophic phase is highly dependent on environmental conditions and is difficult to reproduce in laboratory conditions. Moreover, because the zig-zag model has been so powerful in understanding the biotrophic phase of type III-dependent phytobacteria, their necrotrophic phase is often underestimated.To achieve a better comprehension of plantbacterial interactions, the hemibiotrophicnature of these interactions should be recognized.


S Corre

Coupled Lattice Boltzmann Method for Numerical Simulations of Fully Coupled Heart and Torso Bidomain System in Electrocardiology

Keywords: Coupled Lattice Boltzmann, Heart-torso coupling, coupled nonlinear parabolic/elliptic PDEs, bidomain models, Ionic models, cardiac electrophysiology

In this work, a modied coupling Lattice Boltzmann Model (LBM) in simulation of cardiac electrophysiology is developed in order to capture the detailed activities of macro-to micro-scale transport processes. The propagation of electrical activity in the human heart through torso is mathematically modeled by bidomain type systems. As transmembrane potential evolves, we take into account domain anisotropical properties using intracellular and extracellular conductivity, such as in a pacemaker or an electrocardiogram, in both parallel and perpendicular directions to the bers. The bidomain system represents multi-scale, sti and strongly nonlinear coupled reaction-diusion models that consist of a set of ordinary dierential equations coupled with a set of partial dierential equations. Due to dynamic and geometry complexity, numerical simulation and implementation of bidomain type systems are extremely challenging conceptual and computational problems but are very important in many real-life and biomedical applications. This paper suggests a modied LBM scheme, reliable, ecient, stable and easy to implement in the context of such bidomain systems. Numerical tests to conrm eectiveness and accuracy of our approach are provided and the propagation of electrophysiological waves in the heart is analyzed.

Introduction

Computational cardiac electrophysiological modeling is now an important eld in applied mathematics. Indeed, nowadays, heart and cardiovascular diseases are still the leading cause of death and disability all over the world. That is why we need to improve our knowledge about heart behavior, and more particularly about its electrical behavior. Consequently we want strong methods to compute electrical uctuations in the torso region to prevent cardiac disorders (as arrhythmias), or to study interactions between brain and heart.

In this work, coupling Lattice Boltzmann Models (LBM) in simulation of cardiac electrophysiology are developed in order to capture the detailed activities of macro-to micro-scale transport processes. The propagation of electrical activity in the human heart is mathematically modeled by bidomain type systems in heart, derived from Ohm's law, and is coupled with a general Laplace equation in surrounding tissues across the torso. As heart is embedded in the torso, the dynamic of electrical potential in torso have to be treated in detail to obtain realistic behavior in myocardium (the muscle tissue of the heart) for several applications such as in a pacemaker or an electrocardiogram. The muscle tissue of the ventricles and atria are composed of bers which is oblong like cylinders and formed by thick and thin protein laments. Each ber is connected to neighboring bers with a step-like surface. Because of these connections, the electrical current passes from cell to cell and its direction is therefore determined by the orientation of the bers. This makes the myocardial tissue anisotropic, and we use intracellular and extracellular conductivity, with two tensors respectively called K ic and K ec , in both parallel and perpendicular directions to those bers to taking into account the anisotropy. So, even if this tissue has a granular nature and is composed of multiple individual cells, we can consider this tissue as continuous anisotropic domain which can be likened to a functional syncytium. This biophysical model of electrical cardiac activity links electrophysiological cell models, at small scales, and myocardial tissue mechanics, metabolism and blood ow at large scales. Indeed, myocardium bers have a higher conductivity than normal muscle bers. That is because the heart cells are surrounded by a plasma membrane which induces a dierent chemical and electrical intracellular and extracellular behavior. There is notably a higher intracellular potassium concentration than in the extracellular uid. This dierence induces an outow of potassium which can be traduced by a dierence of distribution charge, called the transmembrane potential ρ. External electrical stimulus control ion channels gates in the membrane, and intracellular and extracellular ionical concentration changes while the transmembrane potential leaves its resting state (see e.g. [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF], [START_REF] Lines | Modeling the electrical activity of the heart: A Bidomain Model of the ventricles embedded in a torso[END_REF] and the references therein).

So, the bidomain system represents multi-scale, sti and strongly nonlinear coupled reaction-diusion models (governing the dynamics at cellular and tissue levels) that consists of a set of ordinary dierential equations coupled with a set of partial dierential equations in myocard to compute extracellular potential ρ e , transmembrane potential ρ and ionic variable U , which are linked to the Laplace equation in torso to compute torso potential ρ T with Robin type boundary conditions between torso and extra-cellular potentials. During last years a lot of studies about bidomain models have led to results about well-posedness, existence and uniqueness of solutions (see e.g., [START_REF] Belmiloudi | Robust control problem of uncertain bidomain models in cardiac electrophisiology[END_REF] , [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF] and the references therein), and several numerical methods based on methods as nite dierence method or nite element method are used to solve these models (see e.g., [START_REF] Sharomi | Convergence order vs. parallelism in the numerical simulation of the bidomain equations[END_REF] and the references therein).

However, there exist several models to simulate electrical behavior in electrophysiology based on Hodgkin-Huxley equations to approximate the electrical characteristics of excitable cells such as neurons and cardiac myocytes. In order to have a simple and realistic model, we choose 2-variable phenomenological Mitchell-Schaeer (MS) model [START_REF] Mitchell | A two-current model for the dinamics of cardiac membrane[END_REF] to capture ionic phenomena. This model derived from Fenton-Karma 3-variable model [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with ber rotation: Filament instability and brillation[END_REF] is interesting because of its simplicity. In the same time, it has already shows its accuracy and its realism. Due to dynamic and geometry complexity, numerical simulation and implementation of bidomain type systems are extremely challenging conceptual and computational problems but are very important in many real-life and biomedical applications. This work presents a modied coupled LBM scheme, reliable, ecient, stable and easy to implement in the context of such bidomain systems with Neumann type boundary conditions in complex geometry boundaries. LBM was originated from Boltzmann's kinetic theory of gases (70s), and attracts more and more attentions for simulating complex uid ows since 90s. More recently, LBM has been extended successfully to simulate dierent types of parabolic reaction-diusion equation as Keller-Segel chemotaxis model [START_REF] Yang | Coupled lattice Boltzmann method for generalized Keller-Segel chemotaxis model[END_REF] and monodomain model in cardiac electrophysiology [START_REF] Campos | Lattice Boltzmann method for parallel simulations of cardiac electrophysiology using GPUs[END_REF], or Poisson equation [START_REF] Chai | Novel Boltzmann model for the Poisson equation[END_REF]. The concept of LBM, which is ecient in terms of parallelization, is based on Boltzmann equations which describe the evolution of particles in kinetic theory. LBM has two main phases: one, local, which models the collision between particles and the second, along each direction of interpolation, which models the transport phase. We demonstrate that solving Lattice Boltzmann Equations (LBE) is equivalent to solving governing bidomain system with Chapman-Enskog perturbation expansion [START_REF] Chapman | The mathematical theory of non-uniform gases[END_REF]. This ansatz leads us to decompose the particular distributions f around equilibrium state in dierent scales of perturbation with introduction of dimensionless Knudsen number ε which is a ratio between the particular mean free path length to a representative physical length scale. So, for every distribution function we can rewrite distribution function f as f = f eq +εf (1) +ε 2 f (2) +.... Then, the main point of Chapman-Enskog expansion is to match LBE with Taylor series with including the decomposition in order to recompose the governing bidomain system. Finally, we propose two applications in simplied square domains for heart and torso. On the one hand, numerical results demonstrate the eectiveness and accuracy of our approach using general methods for bidomain type systems with full known data and chosen analytical solution. On the other one, we treat the model with realistic data to show good agreement with numerical results and potential behavior reported in the literature (see e.g., [START_REF] Boulakia | Mathematical Modeling of Electrocardiograms: A Numerical Study[END_REF], [START_REF] Rioux | A Predictive Method Allowing the use of a Single Ionic Model in Numerical Cardiac Electrophysiology[END_REF]).

The paper is structured as follows: rst, in Section 2 we introduce the full coupled heart and torso bidomain system and the MS ionic model. Then, in Section 3 after some generalities about LBM, we develop our modied coupled LBM for numerical simulations of the considered model and investigate its asymptotic behavior as goes to zero. The validity of this method is demonstrated in Section 4 by comparing the numerical solution to the known exact solution, and convergence of solution is established. Finally, the propagation of electrophysiological waves in the heart is studied. The full coupled heart and torso bidomain model of cardiac tissue is expressed mathematically by the following time dependent system of coupled partial dierential equations governing the electrical potentials in the physical region Ω = Ω H ∪ Ω T , which is an open, bounded, and connected subset of R d , d ≤ 3 and during a time interval (0, T ). The two spaces Ω H and Ω T are the regions occupied by the excitable myocardium tissue (the heart muscle) and the thorax, respectively.

The bidomain Model

First, we will introduce the bidomain model in myocard Ω H . Propagation of electrical potential is modeled by:

div(K ic ∇ρ i ) = I m -κf is , div(K ec ∇ρ e ) = -I m -κf es , (1)
where ρ i and ρ e are the intracellular and extracellular potentials, respectively; K ic and K ec are the conductivity tensors describing the anisotropic intracellular and extracellular conductive media; f is (x, t) and f es (x, t) are the respective externally applied current sources after they passed through the thorax. The transmembrane current density is described by I m and is given by the following expression:

I m = κ(c m ∂ρ ∂t + I ion ), (2)
where ρ is the transmembrane potential, which is dened as ρ = ρ i -ρ e , κ is the ratio of the membrane surface area to the volume occupied by the tissue, c m term is the transmembrane capacitance time unit area. We can assimilate this membrane to a simple electrical circuit with a resistor associate to the ionic current (I ion ) and a capacitor associate to the capacitive current (I c = c m ∂ρ ∂t ) in parallel (Figure 1). The tissue is assumed to be passive, so the capacitance c m can be assumed to be not a function of the state variables. The nonlinear operator I ion (x, t; ρ, U ) describes the sum of transmembrane ionic currents across the cell membrane with U the electrophysiological ionic state variable (which describes e.g., the dynamics of ion-channel and ion concentrations in dierent cellular compartments). These variables satisfy the following ODE (with H ion a nonlinear operator)

∂U ∂t = H ion (x, t; ρ, U ). (3)
From ( 1), ( 2) and ( 3), the bidomain model can be formulated in terms of the state variables ρ, ρ e and U as follows (in

Q H = Ω H × (0, T )) κ(c m ∂ρ ∂t + I ion (.; ρ, U )) -div(K ic ∇ρ) = div(K ic ∇ρ e ) + κf is , -div((K ec + K ic )∇ρ e ) = div(K ic ∇ρ) + κ(f es + f is ), ∂U ∂t = H ion (.; ρ, U ). (4)
The operators I ion and H ion which describe electrophysiological behavior depends on the considered cell ionic model. We can choose physiological ionic models (e.g Djabella and Sorine [START_REF] Djabella | A two-variable model of cardiac action potential with controlled pacemaker activity and ionic current interpretation[END_REF] or Luo and Rudy [START_REF] Luo | A dynamic model of the cardiac ventricular action potential: I. simulations of ionic currents and concentration changes[END_REF]), phenomenological models (e.g Fenton-Karma [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with ber rotation: Filament instability and brillation[END_REF], Fitzhugh-Nagumo [START_REF] Fitzhugh | Impulses and physiological states in theorical models of nerve membrane[END_REF], or Mitchell-Schaeer [START_REF] Mitchell | A two-current model for the dinamics of cardiac membrane[END_REF]). In our case, we choose the phenomenological two-variable classical model proposed by Mitchell and Schaeer as:

I ion (x, t; ρ, U ) = I 1 (x, t; ρ) + I 2 (x, t; ρ)U, H ion (x, t; ρ, U ) = λ 1 (x)U + H 0 (x, t) if ρ ≤ ρ gate , λ 2 (x)U if ρ > ρ gate , (5)
where operators I 1 , I 2 , functions λ 1 , λ 2 and H 0 depend of given parameters and where ρ gate is a given activation threshold.

Then, under the assumption that the heart is isolated from its surrounding (e.g. noncoupled with the thorax), we will search ρ T which is electrical potential in thorax with solving propagation equation (in

Q T = Ω T × (0, T )): -div(K T ∇ρ T ) = f T , (6)
where K T is the conductivity tensor describing the anisotropic conductivity in thorax, and f T is the applied current source.

To close the system, we impose the following boundary conditions on

Σ = ∂Ω H ×(0, T ) (K ic ∇(ρ + ρ e )).n = ζ i , (K ec ∇ρ e ).n = ζ e , (7)
and on Σ T = Γ T × (0, T ) where Γ T = ∂Ω T \ ∂Ω H :

(K T ∇ρ T ).n = 0, (8)
where n being the outward normal to Γ = ∂Ω H and ζ i and ζ e are the intra-and extra-cellular currents per unit area applied across the boundary, and the following initial conditions (in Ω H )

ρ(t = 0) = ρ 0 , U (t = 0) = U 0 .
(9) Such problems have compatibility conditions determining whether there are any solutions to the PDEs. This is easily found by integrating the second equation of (4) over the domain and using the divergence theorem with the boundary conditions (7) (a.e. in (0, T)). Then (for compatibility reasons), we require the following condition

Γ (ζ i + ζ e )dΓ + κ Ω H (f es + f is )dx = 0. (10)
Moreover, the function ρ e is dened within a class of equivalence, regardless of a timedependent function. This function can be xed, for example by setting the following condition (a.e. in (0, T))

Ω H ρ e dx = 0. (11)
Under some hypotheses for the data and parameters of the system and some regularity of operators I ion and H ion , system (4) with ( 7)-( 9) and under the conditions (10)-( 11) is a well-posed problem (for more details see [START_REF] Belmiloudi | Robust control problem of uncertain bidomain models in cardiac electrophisiology[END_REF][START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF]).

We respect the following boundary conditions on Σ:

(K ic ∇ρ + (K ic + K ec )∇ρ e ).n = (K T ∇(ρ T )).n, (12)
which induce the following boundary condition:

(K T ∇ρ T ).n = ζ e + ζ i on Σ. ( 13
)
Remark 2.1 .

In order to guarantee the well-posedness of system (4)-( 9) (under some hypotheses for data and some regularity of the nonlinear operators), we can use the following regularized version of ionic operators H ion , given by ( 5), (see e.g., [START_REF] Boulakia | A Coupled System of PDEs and ODEs arising in Electrocardiograms Modelling[END_REF])

H r (x, t; ρ, U ) = (-λ 2 (x) + (λ 2 (x) -λ 1 (x))h ∞ (ρ)) (U (x, t) -h ∞ (ρ)) -(H 0 (x, t) + λ 1 (x))h ∞ (ρ),
where h ∞ (ρ) is given by

h ∞ (ρ) = 1 2 1 -tanh ρ -ρ gate r g ,
with r g a positive parameter. We can deduce that

lim rg→0 h ∞ (ρ) = 1 if ρ < ρ gate , 0 if ρ > ρ gate ,
To illustrate the convergence of H r to H ion , we consider, e.g., the following functions

λ 1 (x) = 1 T open , λ 2 (x) = 1 T close , H 0 (x, t) = - 1 T open (ρ max -ρ min ) 2 .
with ρ gate = -67mV , T open = 100ms, T close = 120ms, and we choose U (x, t) = 10 -3 . As Figure 2 shows for dierent parameters r g , lim rg→0 H r (ρ, U ) = H ion (ρ, U ).

Numerical method and algorithm

In this section, a numerical method is presented for the bidomain system (4) in two space dimensions. For this, we introduce a coupled modied LBM for solving the coupled system of nonlinear parabolic and elliptic equations (i.e. the rst and second equations of ( 4)). Then we treat the ODE satised by ionic state by applying Gronwall Lemma to obtain an integral formulation, and by using a quadrature rule to approximate the obtained integral.

In the sequel, without loss of generality, we assume c m = 1 and κ = 1. Moreover we assume

K ic = K ic I d , K ec = K ec I d , K T = K T I d , with I d identity matrix and K ic , K ec , K T constants.

LBM for coupled parabolic and elliptic equations

In this rst part, we develop and describe the modied LBM to solve the following system (which corresponds to two rst parts of (4) and elliptic equation in thorax ( 6))

∂ρ ∂t

-div(K ic ∇(ρ + ρ e )) = F (.; ρ, ρ e ), in Q H -div(K ic ∇ρ + (K ic + K ec )∇ρ e ) = G(.; ρ, ρ e ), in Q H -div(K T ∇ρ T ) = H(.; ρ T ), in Q T (14)
where F , G and H are given nonlinear operators.

Before introducing the coupled LBM for system [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF], it is important to explain the main steps of this method and to give some denitions about discretization and operators in case of a general reaction-diusion equation.

LBM for general reaction-diusion equations

Let us introduce the LBM for solving the following reaction-diusion equation (with the macroscopic variable Θ) ∂Θ(x, t) ∂t -div(K(x, t)∇Θ(x, t)) = Φ(x, t; Θ), [START_REF] Corre | Coupled Lattice Boltzmann Modeling of Bidomain type models in Cardiac Electrophysiology[END_REF] where Φ is a non linear operator and K is a diusion coecient. The evolution equation of the LBM for [START_REF] Corre | Coupled Lattice Boltzmann Modeling of Bidomain type models in Cardiac Electrophysiology[END_REF] is given by the following continuous Boltzmann equation [START_REF] Boltzmann | Weitere Studien uber das Warme gleichgenicht unfer Gasmolakuler[END_REF]:

∂φ(x, t; e) ∂t + e • ∇φ(x, t; e) = J (x, t; e) [START_REF] Djabella | A two-variable model of cardiac action potential with controlled pacemaker activity and ionic current interpretation[END_REF] with

Θ(x, t) = φ(x, t; e)de, (17)
where J (x, t; e) = Q col (φ(x, t; e)) + P (x, t; e), φ(x, t; e) is the distribution function of the single particle moving with velocity e at position x and time t, P is the distribution type function of particle of macroscopic external force Φ moving with velocity e and Q col is the Bhatnagar-Gross-Krook (BGK) collision operator dened by (see [START_REF] Bhatnagar | A model for Collision Processes in Gases Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]):

Q col (φ(x, t; e)) = - 1 τ (x, t) (φ(x, t; e) -φ eq (x, t; e)) , (18)
where φ eq is the Maxwell-Boltzmann equilibrium distribution function and τ is the dimensionless relaxation time. LBM leads us to approximate [START_REF] Djabella | A two-variable model of cardiac action potential with controlled pacemaker activity and ionic current interpretation[END_REF] to recover reaction-diusion equation ( 15) with Chapman-Enskog expansion. In the method of Chapman-Enskog the distribution function φ is expanded as : φ = φ (0) + φ (1) + 2 φ (2) + .... = ∞ k=0 k φ (k) , which can be regarded as a power series in a small variable or an expression that keeps track of the relative orders of magnitude of the dierent terms through the scaling parameter (the so-called Kundsen number). The distribution functions φ (0) , φ (1) and φ (2) represent the zero, rst and second approximation to the distribution function φ, and so on. It is important to note that φ eq is an exact solution of the Chapman-Enskog's zero-order approximation of the Boltzmann equation (see [START_REF] Huang | Statistical Mechanics[END_REF]).

The numerical solution of Boltzmann equation ( 16) requires to discretize the conguration spaces, velocities and time. For that, we discretize Q = Ω × (0, T ) in time and space. Then we consider a regular lattice L with typical mesh size ∆x and time step size ∆t, and we dene streaming lattice speed c = ∆x/∆t and lattice sound speed C s by 3C 2 s = c 2 . We note (e i ) i=0,N -1 the discrete set of N admissible particle velocities (the magnitude of each e i depends on the speed c) and we assume that for each node x of lattice L, and each e i , the point x + e i ∆t is also a node of the lattice L. For the latter, we use the so-called D2QN LBM scheme (i.e., 2-dimensional and N velocity vectors). Remark 3.2 .

Since the velocities e i are constant vectors in space-time variables, then, for suciently regular scalar function Ξ, we have

e i • ∇Ξ = div (e i Ξ) . ( 19
)
First, we consider the nite discrete-velocity system of the Boltzmann equation with the nite discrete velocity e i , for i = 0, N -1 (by discretizing the velocity space by the discrete set of microscopic velocities)

∂φ i (x, t) ∂t + e i • ∇φ i (x, t) = J i (x, t) (20)
where, for each particle on the lattice, we associate discrete functions, for i = 0, N -1,

φ i (x, t) = φ(x, t; e i ), φ eq i (x, t) = φ eq (x,
t; e i ) and J i (x, t) = J (x, t; e i ), and the discrete operator Φ i of Φ, which describes the probability of streaming in one particular direction.

The aim of the so-called LBM scheme is to compute the distribution function of particle (φ i (x, t)) i=1,N for x ∈ L and discrete values of time t by solving the discretization of the Boltzmann's equation [START_REF] Zheng | Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method[END_REF] in two steps (for each time): collision and streaming processes. Then, by using a second order time integration scheme to approximate [START_REF] Zheng | Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method[END_REF] (by limiting physical space to a lattice and velocity space to the discrete set of microscopic velocities), we can obtain

φ i (x + e i ∆t, t + ∆t) = φ i (x, t) - 1 τ (x, t) (φ i (x, t) -φ eq i (x, t)) +∆tΦ i (x, t) + ∆t 2 2 ( ∂ ∂t + e i • ∇)Φ i (x, t), (21)
which can be decomposed into the two key sub-steps i.e. collision and streaming. The collision process is local in space and is given for each node x by:

φ col i (x, t) = φ i (x, t)- 1 τ (x, t) (φ i (x, t)-φ eq i (x, t)) +∆tΦ i (x, t)+ ∆t 2 2 ( ∂ ∂t + e i • ∇)Φ i (x, t) (22)
Then we can formally summarize the streaming step as a translation in post-collision distribution function elds. For each node, the post-collision distribution function φ col i (x, t)

will replace post-collision distribution function of the neighbor node φ col i (x + e i ∆t, t) according to velocity e i (Fig. 4). So, streaming process is given for each velocity e i and in a step time ∆t by: φ i (x + e i ∆t, t + ∆t) = φ col i (x, t) (the propagation of Lattice Boltzmann scheme). ( 23) Finally, we recover governing macroscopic variable Θ by a summation of microscopic particle distribution function φ i , i = 1, N as follows (which is a discrete approximation of the integral given in ( 17)) N i=0 φ i (x, t) = Θ(x, t), for x ∈ L and discrete values of time t. [START_REF] Lines | Modeling the electrical activity of the heart: A Bidomain Model of the ventricles embedded in a torso[END_REF] In this paper, we consider the popular two-dimensional D2Q9 lattice, which involves 9 velocity vectors (i.e. N = 9). The microscopic velocities e i , for i = 0, 8 in D2Q9 are given as (see Figure 3)

e 0 = 0 0 , e i = c cos (i -1) π 2 sin (i -1) π 2 , for i = 1, 2, 3, 4, e i = √ 2c cos (i -9 2 ) π 2 sin (i -9 2 ) π 2 , for i = 5, 6, 7, 8.
According to the expression of the microscopic velocities e i , we prove easily that

8 i=0 e i = 8 i=0 w i e i = 0 0 (25)
and

8 i=0 w i e i e i = C 2 s I d , (26)
where the weight vector w = (w 0 , . . . , w 8 ) is dened by (

8 i=0 w i = 1, w 0 = 4/9,
) 27
Remark 3.3 .

With using D2Q9 model, we dene a 9-speed square lattice with square elements. It is also possible to use D2Q7 model to obtain a 7-speed hexagonal lattice with triangular elements or to dene other types of quadrature if it is needed. In all cases, we have to care about geometric properties. Indeed, respect summations [START_REF] Luo | A dynamic model of the cardiac ventricular action potential: I. simulations of ionic currents and concentration changes[END_REF] and ( 26) is necessary to obtain a convergent method. During streaming and collision processes, in order to satisfy boundary conditions, the boundary nodes need special treatments on distribution functions, which are essential to stability and accuracy of the method. We use extrapolation methods [START_REF] Zheng | Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method[END_REF] and modied bounce-back methods [START_REF] Chai | General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method[END_REF], [START_REF] Ladd | Lattice-Boltzmann Simulations of Particle-Fluid Suspensions[END_REF], [START_REF] Yin | An improved bounce-back scheme for complex boundary conditions in lattice Boltzmann method[END_REF].

LBM for coupled system

In this part, we detail the modied coupled LBM presented in [START_REF] Corre | Coupled Lattice Boltzmann Modeling of Bidomain type models in Cardiac Electrophysiology[END_REF], adapted from [START_REF] Yang | Coupled lattice Boltzmann method for generalized Keller-Segel chemotaxis model[END_REF] with introducing LBM for elliptic equations proposed in [START_REF] Chai | Novel Boltzmann model for the Poisson equation[END_REF]. We begin with the introduction of three Lattice Boltzmann type equations, which are associated to the three equations of the coupled system of parabolic and elliptic equations ( 14), respectively. Then, we dene discrete equilibrium functions, source terms, nonlinear operators and relaxation time associated to each equation as in previous section. So, we associate the repartition functions of particle moving with velocity e i , with i = 0, . . . , 8, f i , g i and h i for ρ, ρ e and ρ T , respectively. The discrete LBM equation for the rst equation of ( 14), which is parabolic, can be given by (at node x and time variable t)

f i (x + e i ∆t, t + ∆t) = f i (x, t) - 1 τ 1 (x, t) (f i (x, t) -f eq i (x, t)) +∆tS 1 i (x, t; ∇ρ, ∇ρ e ) +∆tF i (x, t) + ∆t 2 2 ( ∂ ∂t + e i • ∇)F i (x, t), (28)
where τ 1 (x, t) is the dimensionless relaxation factor, F i (x, t) is the discrete source term associated to operator F (x, t; ρ, ρ e ), S 1 i (x, t; ∇ρ, ∇ρ e ) is a link function (between the rst and second equations of ( 14)) related to gradients ∇ρ and ∇ρ e .

For the two elliptic equations of ( 14), for a xed time t, the discrete pseudo-timestepping LBM equations can be given (at node x and pseudo time variables r and s)

g i (x + e i ∆r, r + ∆r; t) = g i (x, r; t) - 1 τ 2 (x, t) (g i (x, r; t) -g eq i (x, r; t)) +∆rS 2 i (x, r; t, ∇ρ, ∇ρ e ) +∆rG i (x, r; t) + ∆r 2 2 ( ∂ ∂r
+ e i • ∇)G i (x, r; t), ( 29
)
h i (x + e i ∆s, s + ∆s; t) = h i (x, s; t) - 1 τ 3 (x, t) (h i (x, t) -h eq i (x, s; t)) +∆sH i (x, s; t) + ∆s 2 2 ( ∂ ∂s
+ e i • ∇)H i (x, s; t), (30)
where ∆r and ∆s are given time-steps, τ 2 (x, t) and τ 3 (x, t) are the dimensionless relaxation factors, g i (x, r; t) and h i (x, s; t) are the distribution functions of particle moving with velocity e i at node x and times r and s, G i (x, r; t) and H i (x, s; t) are the discrete source terms associated respectively to operators G(x, r; t, ρ, ρe ) and H(x, s; t, ρT ). The term S 2 i (x, r; t, ∇ρ, ∇ρ e ) is the second link function (between the rst and second equations of ( 14)) related to gradients ∇ρ and ∇ρ e .

Then, after each time step, we obtain new particular repartitions which lead us to recover governing macroscopic potentials ρ, ρe and ρT with following summations:

8 i=0 f i (x, t) = ρ(x, t), ( 31
) 1 1 -w 0 8 i=0 g i (x, r; t) = ρe (x, r; t), ( 32
) 1 1 -w 0 8 i=0 h i (x, s; t) = ρT (x, s; t). ( 33
)
The two elliptic equations of system (14) associated to the pseudo time variables r and s induce subroutines in order to reach the following convergence:

lim r→∞ ρe (x, r; t) = ρ e (x, t), ( 34
)
lim s→∞ ρT (x, s; t) = ρ T (x, t), (35)
lim r→∞ G(x, r; t, ρ, ρe ) = G(x, t; ρ, ρ e ) [START_REF] Yin | An improved bounce-back scheme for complex boundary conditions in lattice Boltzmann method[END_REF] and lim s→∞ H(x, s; t, ρT ) = H(x, t; ρ T ). [START_REF] Chai | General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method[END_REF] Now we dene equilibrium functions, source terms functions, linked functions and the three dimensionless factors τ 1 , τ 2 , τ 3 . The local equilibrium distribution functions are given by (∀i = 0, . . . , 8)

f eq i (x, t) = w i ρ(x, t), (38)
g eq i (x, r; t) = w i ρe (x, r; t) for i = 0, (w 0 -1)ρ e (x, r; t) for i = 0 (39)

h eq i (x, s; t) = w i ρT (x, s; t) for i = 0, (w 0 -1)ρ T (x, s; t) for i = 0. ( 40
)
According to ( 25) and ( 26), we can prove easily that

8 i=0 f eq i (x, t) = ρ(x, t), (41)
8 i=0 e i f eq i (x, t) = 0 0 , ( 42
) 8 i=0 e i e i f eq i (x, t) = C 2 s ρ(x, t)I d , (43)
8 i=0 g eq i (x, r; t) = 0, (44) 8 i=0
e i g eq i (x, r; t) = 0 0 ,

e i e i g eq i (x, r; t) = C 2 s ρe (x, r; t)I d , (46)

8 i=0 h eq i (x, s; t) = 0, (47)
8 i=0 e i h eq i (x, s; t) = 0 0 , (48) 8 i=0
e i e i h eq i (x, s; t) = C 2 s ρT (x, s; t)I d .

(49) Remark 3.5 .

There is no uniqueness of equilibrium function denition. According to implementation choices and type of PDE, we can choose more complex forms which are more ecient in term of implementation. This choice only depends on the denitions of link functions

S 1
i and S 2 i , and time relaxations τ 1 and τ 2 . Indeed, we can take into account nonlinear diusion operators in relaxation terms τ 1 and τ 2 , and we can merge equilibrium and link functions (see e.g. [START_REF] Xiang | modied lattice boltzman scheme for nonlinear convection diusion equations[END_REF]).

Then, we can dene the source terms F i , G i and H i , which are discrete approximations of nonlinear operator and source terms of system [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF], as follows (∀i = 0, . . . , 8)

F i (x, t) = w i F (x, t; ρ, ρ e ), (50)
G i (x, r; t) = w i G(x, r; t, ρ, ρe ), (51)
H i (x, s; t) = w i H(x, s; t, ρT ).
(52)

According again to [START_REF] Luo | A dynamic model of the cardiac ventricular action potential: I. simulations of ionic currents and concentration changes[END_REF] and [START_REF] Mitchell | A two-current model for the dinamics of cardiac membrane[END_REF], we can derive the following relations

8 i=0 F i (x, t) = F (x, t; ρ, ρ e ), ( 53
) 8 i=0 e i F i (x, t) = 0 0 , (54) 8 i=0
G i (x, r; t) = G(x, r; t, ρ, ρe ),

(55)

8 i=0 e i G i (x, r; t) = 0 0 , ( 56
) 8 i=0 H i (x, t) = H(x, s; t, ρT ), (57) 8 i=0
e i H i (x, s; t) = 0 0 .

(58)

Similarly, we dene the rst link function S 1 i associate to the rst equation of ( 14) and equation ( 28) by (∀i = 0, . . . , 8)

S 1 i (x, t; ∇ρ, ∇ρ e ) = w i (β ρ (x, t)e i • ∇ρ e (x, t) + α ρ (x, t)e i • ∇ρ(x, t)), (59)
and the second link function S 2 i associate to the second equation of ( 14) and equation ( 29) as (∀i = 0, . . . , 8)

S 2 i (x, r; t, ∇ρ, ∇ρ e ) = w i (β ρe (x, t)e i • ∇ρ(x, t) + α ρe (x, t)e i • ∇ρ e (x, r; t)), (60)
where we choose α ρ , α ρe , β ρ , β ρe to be:

α ρ (x, t) = - 1 τ 1 (x, t)C 2 s ∆t K ic (x, t) + C 2 s ∆t 1 2 -τ 1 (x, t) , (61)
β ρ (x, t) = -K ic (x, t) τ 1 (x, t)C 2 s ∆t , (62)
α ρe (x, t) = - 1 τ 2 (x, t)C 2 s ∆r (K ic (x, t) + K ec (x, t)) + C 2 s ∆r 1 2 -τ 2 (x, t) , (63)
βρ e (x, t) = -K ic (x, t) τ 2 (x, t)C 2 s ∆r . ( 64
)
Remark 3.6 .

As values α ρ and α ρe dened by the equations ( 61) and (64) depend on relaxation factors τ 1 and τ 2 , we can choose them as arbitrary nonzero constants. However, for the relaxation time τ 3 , in order to recover relaxation tensor K T I d in the governing elliptic equation, we can dene it as

τ 3 (x, t) = K T (x, t) C 2 s ∆s + 1 2 .
(65)

By using the relations ( 25) and ( 26), we can deduce the following properties Remark 3.7 .

The link functions S 1 i and S 2 i depend on the gradient of functions ρ, ρ e and ρe . The approximation of these gradients can be directly derived from Chapman-Enskog expansion with distribution functions f i and g i (see the below equations ( 101) and ( 111)), or with a nite dierence approximation type method.

We can now prove formally that the above coupled LBM recovers the governing macroscopic system [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF]. Proposition 3.1 .

If local equilibrium distribution functions f eq i , g eq i and h eq i satisfy relations (41)-(49), source terms F i , G i , H i respect properties given by equations ( 53)-(58), properties (66)-(69) hold for link functions S 1 i and S 2 i , then, for

τ 3 = K T C 2 s ∆t + 1
2 , the above coupled LBM recover correctly macroscopic system [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF] through the Chapman-Enskog analysis.

Proof: To derive the macroscopic equation from developed coupled lattice BGK model, the Chapman-Enskog expansion is applied under the assumption of small Kundsen number (which can be a ratio between a characteristic length L and the particular mean free path l) to determine equilibrium functions, source terms and link functions. Each function is decomposed around its equilibrium state with dierent scale of perturbations. Then, the dierence between associated LBE and Taylor series of the distribution function lead us to recover the governing macroscopic equation.

Before to prove the main result of proposition, we introduce the main steps of demonstration with applying Chapman-Enskog expansion to distribution functions φ i of LBE [START_REF] Huang | A Fully Implicit Method for Lattice Boltzmann Equations[END_REF] in order to recover governing reaction-diusion equation [START_REF] Corre | Coupled Lattice Boltzmann Modeling of Bidomain type models in Cardiac Electrophysiology[END_REF]. The macroscopic variable Θ is dened in terms of distribution functions as e i e i φ eq i (x,

t) = C 2 s Θ(x, t)I d . (70)
For every distribution functions φ i the Chapman-Enskog expansion is applied with the small parameter as:

φ i = φ (0) i + εφ (1) i + ε 2 φ (2) i + ... = φ eq i + εφ (1) i + ε 2 φ (2) i + ..., (71)
∂ ∂t = ε ∂ ∂t 1 + ε 2 ∂ ∂t 2 , (72)
∇ = ε∇ 1 , div 1 = ∇ 1 •, ( 73
)
where t i is the time scale, for i = 1, 2, x 1 is the space scale and ∇ 1 = ∇ x1 is the gradient with respect to x 1 , and we assume that every source term Φ i take the form

Φ i = εΦ (1) i + ε 2 Φ (2)
i , with

8 i=0 Φ (n) i (x, t) = Φ (n) (x, t) and 8 i=0 e i Φ (n) i (x, t) = 0 0 . ( 74
)
Remark 3.8 .

The Chapman-Enskog expansion is usually called multi-scale Chapman-Enskog expansion, because this ansatz lead us to decompose every function according to several perturbation scales. But we take into account the multi-scale nature of a biological or physical problem through the relaxation-time τ (x, t) which is related to diusion properties and nondimensionalization parameters.

Property 3.1 .

The functions φ (0) i are associated to unperturbed state (i.e.: ∀i = 0, . . . , 8 : φ (0) i = φ eq i ). Moreover, since each scale of perturbation is linearly independent of each other, we can deduce that (for all n ∈ IN * )

8 i=0 φ (n) i = 0.
(75) Then, we express the dierence between [START_REF] Huang | A Fully Implicit Method for Lattice Boltzmann Equations[END_REF] and Taylor series of distribution function φ i moving with velocity e i . Indeed Taylor series also lead us to approximate the term φ i (x + e i ∆t, t + ∆t) as:

φ i (x + e i ∆t, t + ∆t) = φ i (x, t) + ∞ k=1 ∆t k k! ∂ ∂t + e i • ∇ k φ i (x, t), (76)
where

∂ ∂t + e i • ∇ k = k n=0 k n ∂ k-n ∂t k-n (e i • ∇) n ,
with

(e i • ∇) n = (e i • ∇) (e i • ∇) n-1 .
So, by using the expressions (74)-(73), Taylor series (76) and discrete LBE ( 21) become (∀i = 0, . . . , 8)

φ i (x + e i ∆t, t + ∆t) = φ (0) i (x, t) + εφ (1) i (x, t) + ... + ∞ k=1 ∆t k k! ε ∂ ∂t 1 + ε 2 ∂ ∂t 2 + e i • (ε∇ 1 ) k φ (0) i (x, t) + εφ (1)
i (x, t) + ...

and

φ i (x + e i ∆t, t + ∆t) = 1 - 1 τ φ (0) i (x, t) + εφ (1)
i (x, t) + ...

+ 1 τ φ eq i (x, t) + ∆t εΦ (1) i (x, t) + ε 2 Φ (2) i (x, t) + ∆t 2 2 ε ∂ ∂t 1 + ε 2 ∂ ∂t 2 + εe i • ∇ 1 εΦ (1) i (x, t) + ε 2 Φ (2) i (x, t) . ( 78
)
By calculating the dierence between relations (78) and (77), we can obtain

0 = 1 τ φ (0) i (x, t) + εφ (1)
i (x, t) + ...

+ ∞ k=1 ∆t k k! ε ∂ ∂t 1 + ε 2 ∂ ∂t 2 + e i • (ε∇ 1 ) k φ (0) i (x, t) + εφ (1)
i (x, t) + ...

- 1 τ φ eq i (x, t) -∆t εΦ (1) i (x, t) + ε 2 Φ (2) i (x, t) - ∆t 2 2 ε ∂ ∂t 1 + ε 2 ∂ ∂t 2 + εe i • ∇ 1 εΦ (1) i (x, t) + ε 2 Φ (2) i (x, t) , (79)
which can be rewritten as follows (according to dierent scales of ε)

0 = Ψ (0) i (x, t) + εΨ (1) i (x, t) + ε 2 Ψ (2) i (x, t) + . . . , (80)
where probability density functions Ψ (n) i regroup all terms of (79) in order of ε n , n = 0, ..., respectively (which are partial dierential equations for n = 0)

Ψ (0) i = 1 τ φ (0) i -φ eq i = 0, (81)
Ψ (1) i = ∆t ∂ ∂t 1 + e i • ∇ 1 φ (0) i + 1 τ φ (1)
i -∆tΦ

(1) i = 0, (82)
Ψ (2) i = ∆t 2 2 ∂ ∂t 1 + e i • ∇ 1 2 + ∆t ∂ ∂t 2 φ (0) i + ∆t ∂ ∂t 1 + e i • ∇ 1 φ (1) i + 1 τ φ (2) i -∆tΦ (2)
i -

∆t 2 2 ( ∂ ∂t 1 + e i • ∇ 1 )Φ (1) i = 0, (83)
. . .

As functions Ψ i are density probability functions, summation along each velocity e i lead us to recover macroscopic functions, noted Ψ (n) at order n. Then, we recover macroscopic functions with

8 i=0 Ψ (0) i + εΨ (1) i + ... = Ψ (0) + εΨ (1) + . . . ( 84
)
and with an order N accuracy as:

Ψ (0) + εΨ (1) + ... + ε N -1 Ψ (N -1) + O(ε N ) = Ψ + O(ε N ). (85)
Finally, by retaining terms up to O(ε 3 ), taking (82) times plus (83) times 2 , summing over i and using (81), ( 70), ( 72), ( 73) and (75), we can obtain, by a suitable choice of equilibrium distribution φ eq i (with similar properties as f eq i , g eq i and h eq i given in ( 41)-(49)), source terms Φ i (with similar properties as F i , G i and H i given in (53)-( 58)) and relaxation time τ (x, t) = K(x,t) C 2 s ∆t + 1 2 , the approximation of the general reaction-diusion equation [START_REF] Corre | Coupled Lattice Boltzmann Modeling of Bidomain type models in Cardiac Electrophysiology[END_REF].

In fact, according to (74),(70), (75) and using the relation [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with ber rotation: Filament instability and brillation[END_REF], we can prove easily from (82) that

Ψ (1) = ∆t( ∂Θ(x, t) ∂t 1 -Φ (1) (x, t)). (86)
To treat the second order terms, we rewrite

∆t 2 2 ∂ ∂t1 + e i • ∇ 1 2
φ eq i with using rst order term given by (82) as:

∆t 2 2 ∂ ∂t 1 + e i • ∇ 1 2 φ eq i = ∆t 2 ∂ ∂t 1 + e i • ∇ 1 - 1 τ φ (1)
i + ∆tΦ

(1) i , = -∆t ∂ ∂t 1 + e i • ∇ 1 1 2τ φ (1) i + ∆t 2 2 ∂ ∂t 1 + e i • ∇ 1 Φ (1) i and then Ψ (2) i becomes Ψ (2) i = ∆t ∂ ∂t 1 + e i • ∇ 1 (1 - 1 2τ )φ (1) i + 1 τ φ (2) i +∆t ∂ ∂t 2 φ eq i -∆tΦ (2)
i .

(87)

Summing over i the terms of equation ( 87), according again to (74),(70), (75) and using the relation [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with ber rotation: Filament instability and brillation[END_REF], we can deduce that

Ψ (2) = ∆t ∂Θ(x, t) ∂t 2 -Φ (2) + div 1 1 - 1 2τ 8 i=0 e i φ (1) i
.

To express 1 -

1 2τ 8 i=0 e i φ (1)
i , we use relation (82) and properties (74),(70) (using ( 19))

1 - 1 2τ 8 i=0 e i φ (1) i = -τ C 2 s ∆t 1 - 1 2τ ∇ 1 Θ.
According to the expression of τ we can deduce that τ C 2 s ∆t 1 -

1 2τ
= K, and then

Ψ (2) = ∆t( ∂Θ(x, t) ∂t 2 -div 1 (K∇ 1 Θ(x, t)) -Φ (2) ). (88)
Then, we can recompose Ψ with using (81), ( 85), ( 86) and ( 88) and deduce the following approximation of equation [START_REF] Corre | Coupled Lattice Boltzmann Modeling of Bidomain type models in Cardiac Electrophysiology[END_REF] as

∂Θ(x, t) ∂t -div(K(x, t)∇Θ(x, t)) -Φ(x, t; Θ) = 0. (89)
Now, we will follow the previous process and steps for each equation of bidomain system ( 14) and its associate discrete LBE system ( 28)- [START_REF] Satofuka | Parallelization of lattice Boltzmann method for incompressible ow computations[END_REF]. For this, we perform the Taylor expansion on ( 28)-( 30) up to terms with order of O(ε 3 ) and we use ansatz given by equations ( 71)-(73) where distribution functions f i , g i and h i successively replace φ i with order 2 expansion given by equation ( 71) and we assume that the source terms and link functions F i ,G i , H i , S 1 i and S 2 i take the form (which are dened like source term Φ i with order 2 expansion)

F i = εF (1) i + ε 2 F (2) i , with 8 i=0 e i F (j) i = 0, for j = 1, 2, (90)
S 1 i = εS 1, (1) i
, with

8 i=0 S 1,(1) i = 0, (91)
G i = εG (1) i + ε 2 G (2) i , with 8
i=0 e i G (j) i = 0, for j = 1, 2, (92)
S 2 i = εS 2, (1) i
, with

8 i=0 S 2,(1) i = 0, (93)
H i = εH (1) i + ε 2 H (2) i , with 8 i=0 e i H (j) i = 0, for j = 1, 2. (94)
First, we consider the rst equation of ( 14) and its associate discrete LBE system [START_REF] Premnath | Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase ow[END_REF].

According to (90)-(91) and using similar arguments to derive the terms 81)-(87), we can obtain the following terms P (n) i , n = 0, 2 (after division by ∆t)

Ψ (n) i in (
P (0) i = - 1 ∆tτ 1 f (0) i -f eq i = 0, (95)
P (1) i = ∂ ∂t 1 + e i • ∇ 1 f eq i + 1 τ 1 ∆t f (1) i -S 1,(1) i -F (1) i = 0, (96)
P (2) i = ∂ ∂t 1 + e i • ∇ 1 f (1) i + ∂ ∂t 2 + ∆t 2 ∂ ∂t 1 + e i • ∇ 1 2 f eq i + 1 τ 1 ∆t f (2) i - ∆t 2 ∂ ∂t 1 + e i • ∇ 1 F (1) i -F (2) i = 0. ( 97
) i=0 f (n) i
= 0, for n = 0, and the relations (41), (42), ( 53), (66), ( 90) and (91) we obtain

P (1) = 8 i=0 P (1) i = ∂ρ ∂t 1 -εF (1) = 0. (98)
By rewriting the equation (97) with the equation (96), we obtain

P (2) i = ∂ ∂t 1 + e i • ∇ 1 (1 - 1 2τ 1 )f (1) i + ∂f eq i ∂t 2 + 1 τ 1 ∆t f (2) i + ∆t 2 ∂ ∂t 1 + e i • ∇ 1 S 1,(1) i -F (2)
i .

(99)

Summing over i the terms of equation ( 99) and using the fact that

8 i=0 f (n) i
= 0, for n = 0, Remark 3.2 and the relations (41), ( 42), ( 53), (66), ( 90) and (91) we obtain

2 P (2) = 2 8 i=0 P (2) i = 2 ∂ρ ∂t 2 + div (1 - 1 2τ 1 ) 8 i=0 e i f (1) i + ∆t 2 div( 8 i=0 e i S 1 i ) -2 F (2) , (100)
with

F (j) = 8 i=0 F (j)
i , for j = 1, 2.

Now we calculate

(1 -1 2τ1 ) 8 i=0 e i f (1)
i . From (96) (to express

8 i=0 e i f (1)
i ) and according to [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with ber rotation: Filament instability and brillation[END_REF], (54), (68), (90) and (91), we can deduce

(1 - 1 2τ 1 ) 8 i=0 e i f (1) i
= -(1 - 1 2τ 1 )C 2 s τ 1 ∆t∇ρ + τ 1 ∆t(1 - 1 2τ 1 ) 8 i=0 e i S 1 i . (101)
According to (68) and (101), the relation (100) becomes

2 P (2) = 2 ∂ρ ∂t 2 -2 F (2) +div C 2 s τ 1 ∆t β ρ ∇ρ e + (α ρ + (1 - 1 2τ 1
))∇ρ .

(102)

Then (from (59), ( 61) and ( 62)

) 2 P (2) = 2 ∂ρ ∂t 2 -2 F (2) -div (K ic (∇ρ + ∇ρ e )) . (103)
Consequently, from (98) and (103), we can recover the macroscopic parabolic equation of ( 14) with

0 = P (0) + εP (1) + ε 2 P (2) + O(ε 3 ) = ∂ρ ∂t -F -div (K ic (∇ρ + ∇ρ e )) . (104)
Second, for a xed time t, we consider the second equation of ( 14) and its associate discrete pseudo-time-stepping LBM system (29) (for pseudo-time variables r with

∂ ∂r = ε ∂ ∂r 1 + ε 2 ∂ ∂r 2
and pseudo-time step ∆r). By using similar arguments to derive the terms

P (n)
i , n = 0, 2, we can obtain the following terms (after division by ∆r)

Q (0) i = - 1 ∆tτ 2 g (0) i -g eq i = 0, (105)
Q (1) i = ∂ ∂r 1 + e i • ∇ 1 g eq i + 1 τ 2 ∆r g (1) i -S 2,(1) i -G (1) i = 0, (106)
Q (2) i = ∂ ∂r 1 + e i • ∇ 1 g (1) i + ∂ ∂r 2 + ∆r 2 ∂ ∂r 1 + e i • ∇ 1 2 g eq i + 1 τ 2 ∆r g (2)
i -

∆r 2 ∂ ∂r 1 + e i • ∇ 1 G (1) i -G (2) i = 0. ( 107
)
Summing over i the terms of equation ( 106) and using the fact that

8 i=0 g (n) i
= 0, for n = 0, and the relations (44), ( 45) and ( 93) we obtain

Q (1) = -G(1) , (108)
with

G(1) = 8 i=0 G (1)
i . By rewriting the equation (107) with the equation (106), we obtain

Q (2) i = ∂ ∂r 1 + e i • ∇ 1 (1 - 1 2τ 2 )g (1) i + ∂g eq i ∂r 2 + 1 τ 2 ∆r g (2) i + ∆r 2 ∂ ∂r 1 + e i • ∇ S 2,(1) i -G (2)
i .

(109)

Summing over i the terms of equation (109) and using the fact that

8 i=0 g (n) i = 0, for n = 0, 8 i=0
g eq i = 0, Remark 3.2 and the relations (44), ( 45)and (93) we obtain

2 Q (2) = 2 8 i=0 Q (2) i = div (1 - 1 2τ 2 ) 8 i=0 e i g (1) i + ∆r 2 div( 8 i=0 e i S 2 i ) -2 G(2) , (110)
with

G(2) = 8 i=0 G (2)
i .

Now we calculate

(1 -1 2τ2 ) 8 i=0 e i g (1)
i . From (106) and according to [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with ber rotation: Filament instability and brillation[END_REF], (54), (69), (92) and (93), we can deduce

(1 - 1 2τ 2 ) 8 i=0 e i g (1) i = -(1 - 1 2τ 2 )C 2 s τ 2 ∆r∇ρ + τ 2 ∆r(1 - 1 2τ 1 ) 8 i=0 e i S 2 i . (111)
According to ( 69) and ( 111), the relation (110) becomes

2 Q (2) = -2 G(2) +div C 2 s τ 2 ∆r β ρe ∇ρ + (α ρe + (1 - 1 2τ 2
))∇ρ e .

(112)

Then (from (60), ( 63) and ( 64))

2 Q (2) = -2 G(2) -div (K ic ∇ρ + (K ic + K ec )∇ρ e ) . ( 113
)
We can now sum the terms of order 0 to 2 to obtain the following equation (from (108) and ( 113))

0 = Q (0) + εQ (1) + ε 2 Q (2) + O(ε 3 ) = -G -div (K ic ∇ρ + (K ic + K ec )∇ρ e ) . ( 114
)
To recover the rst elliptic equation of system [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF], we use the two limits ( 34) and (36) to write the following limit: (formally)

lim r→∞ -G(x, r; t, ρ, ρ e ) -div (K ic ∇ρ(x, t) + (K ic + K ec )∇ρ e (x, r; t)) = -G(x, t; ρ, ρ e ) -div (K ic ∇ρ(x, t) + (K ic + K ec )∇ρ e (x, t)) . ( 115
)
Finally, we apply the previous method to recover the governing second elliptic equation of the system [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF]. For that, for a xed time t, we consider the discrete pseudo-timestepping LBM system (30) (for pseudo-time variables s with

∂ ∂s = ε ∂ ∂s 1 + ε 2 ∂ ∂s 2
and pseudo-time step ∆s). By using similar arguments to derive the terms Q (n) i , n = 0, 2, we can obtain the following terms (after division by ∆s)

R (0) i = - 1 τ 3 ∆s h (0) i -h eq i = 0, (116)
R (1) i = ∂ ∂s 1 + e i • ∇ 1 h eq i + 1 τ 3 ∆s h (1) i -H (1) i = 0, (117)
R (2) i = ∂ ∂s 1 + e i • ∇ 1 h (1) i + ∂ ∂s 2 + ∆s 2 ∂ ∂s 1 + e i • ∇ 1 2 h eq i + 1 τ 3 ∆s h (2) i - ∆s 2 ∂ ∂s 1 + e i • ∇ 1 H (1) i -H (2) i = 0. ( 118
)
Summing over i the terms of equation ( 117) and using the fact that

8 i=0 h (n) i
= 0, for n = 0, and the relations ( 47) and ( 48) we obtain

εR (1) = -ε H(1) , (119)
with

H(1) = 8 i=0 H (1)
i . By rewriting the equation (118) with the equation (117), we obtain

R (2) i = ∂ ∂s 1 + e i • ∇ 1 (1 - 1 2τ 3 )h (1) i + ∂h eq i ∂s 2 + 1 τ 3 ∆r h (2) i -H (2)
i .

(120

) i=0 h (n) i = 0, for n = 0, 8 i=0
h eq i = 0, Remark 3.2 and the relations (47) and ( 48) we obtain

2 R (2) = 2 8 i=0 R (2) i = div (1 - 1 2τ 3 ) 8 i=0 e i h (1) i -2 H(2) , (121)
with

H(2) = 8 i=0 H (2)
i .

Now we calculate

(1 -1 2τ3 ) 8 i=0 e i h (1)
i . From (117) and according to [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with ber rotation: Filament instability and brillation[END_REF], ( 58) and (65), we can deduce

(1 - 1 2τ 3 ) 8 i=0 e i h (1) i = -(1 - 1 2τ 3 )C 2 s τ 3 ∆s∇ρ T = -K T ∇ρ T . ( 122
)
According to (122), the relations (121) becomes

2 R (2) = -2 H(2) -div (K T ∇ρ T ) . ( 123
)
We can now sum the terms of order 0 to 2 to obtain the following equation (from ( 119) and ( 123))

0 = R (0) + εR (1) + ε 2 R (2) + O(ε 3 ) = -H -div (K T ∇ρ T ) . ( 124
)
To recover the second elliptic equation of system ( 14), we use the two limits ( 35) and (37) to write the following limit (formally):

lim s→∞ -H(x, s; t, ρ T ) -div (K T ∇ρ T (x, s; t)) = -H(x, t; ρ T ) -div (K T ∇ρ T (x, t)) . ( 125
)
This completes the proof.

Treatment of ODE

Now, we present briey the method to solve ODE satisfy by variable V on the interval [s 1 , s 2 ], with initial condition V (x, s 1 ) = V 0 (x). As:

∂V (x, s) ∂s = W (x, s; V ), (126)
where operator W is dened by

W (x, s; V ) = W (x, s) -λ(x)V (x, s). ( 127
)
By using Gronwall Lemma we can deduce that

V (x, s 2 ) = V 0 (x)e -λ(x)(s2-s1) + e -λ(x)s1 s2 s1
W (x, s)e λ(x)s ds.

Then, according to approximation of derived integral by trapezoidal method between t 1 and t 2 we obtain the following approximation of V denoted also by V V (., s 2 ) = V 0 (.)e -λ(.)(s2-s1) + s 2 -s 1 2 W (., s 1 ) + W (., s 2 )e λ(.)(s2-s1) .

(128)

In our case, we assume that for all t ∈ [0; T -∆t] the ionic function U given by (3) looks like V given by (127) between t and t + ∆t . So with thin enough discretization, we aim to reach an order 2 convergence space accuracy method to approximate ionic variable U with using relation (128).

3.3

Algorithm

After non-dimensionalization, mesh denition and initialization of initial conditions, parameters and data, the proposed algorithm to solve the full coupled heart and torso bidomain system can be summarized as follows 1. Initialization: t=0. 2. LBE according to time t by using [START_REF] Premnath | Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase ow[END_REF] to compute ρ(x, t + ∆t).

3. Trapezoidal method by using (128) to compute U (x, t + ∆t).

Loop on new time variable r:

(a) LBE according to time r by using [START_REF] Rioux | A Predictive Method Allowing the use of a Single Ionic Model in Numerical Cardiac Electrophysiology[END_REF] to compute ρe (x, r + ∆r; t + ∆t) (i.e., approximation of ρ e (x, t + ∆t)) . (b) If we don't reach the stopping criteria chosen to be:

|g i (x, r + ∆r; t) -g i (x, r; t)| < C ρe ∆x Ne , ( 129
)
where C ρe is a chosen constant and N e the order of accuracy we want, set r = r + ∆r and go back to (4a).

5. Set ρ e (x, t + ∆t) := ρe (x, r; t + ∆t).

Loop on new time variable s:

(a) LBE according to time s by using [START_REF] Satofuka | Parallelization of lattice Boltzmann method for incompressible ow computations[END_REF] to compute ρT (x, s + ∆s; t + ∆t) (i.e., approximation of ρ T (x, t + ∆t)), with taking into account boundary condition ( 12) between ρ T and ρ e . (b) If we don't reach the stopping criteria chosen to be:

|h i (x, s + ∆s; t) -h i (x, s; t)| < C ρ T ∆x N T , ( 130
)
where C ρ T is a chosen constant and N T the order of accuracy we want, set s = s + ∆s and go back to (6a).

7. Set ρ T (x, t + ∆t) := ρT (x, s; t + ∆t).

8. If t = T , set t := t + ∆t and go back to (2). Remark 3.9 .

In order to parallelize the LBM, the main idea is to treat separately the collision and the streaming phases. The collision step is strictly local and can be computed on each node independently of its neighbors. Then, the streaming is just an exchange of data between its node of origin and the adjacent neighbor along velocity (see e.g. [START_REF] Satofuka | Parallelization of lattice Boltzmann method for incompressible ow computations[END_REF], [START_REF] Zhi | Parallelization of Lattice Boltzmann method using MPI domain decomposition technology for a drop impact on a wetted solid wall[END_REF], and references therein).

Numerical simulation and applications

To validate the capacity of our modied coupled LBM to deal with 2D bidomain systems, several situations are numerically simulated. In this study, we consider the phenomenological two-variable model proposed by Mitchell and Schaeer [START_REF] Mitchell | A two-current model for the dinamics of cardiac membrane[END_REF]. The functions H ion and I ion are given by:

I ion = - U T in (ρ -ρ min ) 2 (ρ max -ρ) ρ max -ρ min + 1 T out ρ -ρ min ρ max -ρ min (131)
and

H ion = U Topen - 1 Topen(ρmax-ρmin) 2 if ρ < ρ gate , U T close if ρ ≥ ρ gate . (132)
Then, the ionic variable U depends on the change-over voltage ρ gate , the resting potential ρ min , the maximum potential ρ max , and on the time constants for opening, T open , and closing, T close . The time constants T in and T out are respectively related to the length of depolarization and repolarization phases. These constants are such that

T in < T out < min(T open , T close ).
In the sequel, we take the values of the parameters, which are consistent with the analysis and observation reported in Mitchell and Schaeer [START_REF] Mitchell | A two-current model for the dinamics of cardiac membrane[END_REF], given e.g. in [START_REF] Boulakia | Mathematical Modeling of Electrocardiograms: A Numerical Study[END_REF] (see Table 1). Moreover, we assume that the domain heart-torso is in a square region Ω = Description name value (unit) Cell surface to volume ratio [-L T /2; L T /2] × [-L T /2; L T /2] = Ω H ∪ Ω T , where Ω H = [0; L H ] × [0; L H ] and Ω T = Ω \ Ω H , where L H and L T are respectively heart and torso lengths (Figure 5). Nota Bene: If the exact solution φ sol is known, we can measure the eciency of method with the following L 2 relative error:

Err φ = φ sol -φ L 2 (Ω) φ sol L 2 (Ω) .

Benchmark problem and validation

In this rst analysis, we investigate the accuracy and spatial convergence rate of the proposed modied LBM for which we postulate that the error estimates of the method is of order 2 in space and of order 1 in time (for suciently regular solution). We assume L H = 1, L T = 3 and T = 1, and we take K ic = 1, K ec = 1, ξ i = 0, ξ e = 0, ρ(x, t = 0) = 0 and U (x, t = 0) = cos(π(x + y)). Moreover the source terms F ρ , F ρe and F ρ T are added to parabolic-elliptic system (4) and elliptic equation in thorax [START_REF] Boltzmann | Weitere Studien uber das Warme gleichgenicht unfer Gasmolakuler[END_REF], respectively, to ensure that the function (ρ, ρ e , U, ρ T ) given by (with x = (x, y)) :

ρ(x, t) = tx 2 (x-1) 2 y 2 (y -1) 2 , ρ e (x, t) = t(cos(πx)+cos(πy)), U (x, t) = e t cos(π(x+y))

ρ T (x, t) = K ic ρ(x, y, t) + (K ic + K ec )ρ e (x, y, t)
is the exact solution of system ( 4)-( 6) (with appropriate source term).

To study the convergence, we have constructed a sequence of meshes with decreasing spatial step ∆x between 1/25 and 1/200 and ∆t = ∆x 2 . Figures 6 ,8, 7, 9 show respectively numerical solution of ρ, ρ e , U and ρ T at t = 1 for the step ∆x = 1/100. Since, the chosen stopping criteria for the iterative method to approach ρ e and ρ T involves a constant error (because this criteria is not dened in function of lattice size cell), then we analyze only the numerical estimate of relative error of approximation of ρ and U .

We present on Figures 10 and11 (at t = 0.5 and t = 1) the convergence curves, log(Error) versus ∆x, for ρ and U (numerical values can be found in Table 2). We observe that the slope of error curves for ρ passes approximately from 1.5 to 2 and the slope of error curves for U is approximately equal to 2. This shows that our numerical error estimates agree with the postulated error estimates, and indicate the good performance of our method. Table 2: Relation between relative error and lattice spacing for ρ and U.

-8 -7.5 -7 -

t=0.5 t=1

Figure 10: Error curves with respect to spatial step ∆x for ρ in ΩT at t = 0.5 and t = 1.

-8 -7.5 -7 - Remark 4.1 .

We notice that the two elliptic equations have relative errors which depend on convergence criteria given in ( 129) and (130) with C ρe = C ρ T = 1 and N T = N e = 2. These criteria lead to a 10 -2 order relative error for both solutions of elliptic equations.

Study of an action potential in heart and torso

With this second set of data, we aim to reproduce a realistic electrical wave propagation in heart and torso. For that, we apply an external density current within source term f is in a short time scale (t ∈ [0, T act ]) which represents external applied current (electroshock for example). This applied current involves action potentials (AP) which are generated by the movement of ions into and out of cardiac cells (i.e. biological manifestation of electrical cells excitation) with diusion in torso region. This potential leads us then to study the evolution of the transmembrane potential ρ and ionic variable U . The purpose of the next analysis is to derive the following standard four phases of the action potential (see e.g. [START_REF] Britton | Essential Mathematical Biology[END_REF] and references therein).

• Phase I: Depolarization. In a short time scale (order of T in ), transmembrane potential (TMP) changes from its resting state ρ min to its maximum ρ max when cells are depolarized. This sudden change induces opening of ionic channels, so ionic variable U leaves its equilibrium state to falling down.

• Phase II: Excited phase. In a long time scale (order of T open ), TMP ρ begins to fall slowly. During this phase, ionic channels are open and ionic variable U goes to its minimum.

• Phase III: Repolarization. In a short time scale (order of T out ), TMP goes back towards resting potential of ρ min to prepare the cells for a new cycle of depolarization. During this phase, ionic channels are gradually inactivated. That induces inversion of ionic variable slope.

• Phase IV: Recovery period. In a long time scale (order of T close ), TMP reaches the resting potential ρ min . During this time, ionic variable slowly goes back to its equilibrium state and ionic channels are closed at resting TMP.

To reproduce this behavior, we take the values of the parameters given in Table 1, we choose discretization parameters ∆x = 1/50, ∆t = ∆x 2 and convergence criteria parameters C ρe = C ρ T = 1 and N e = N T = 3, we x L H = 15cm and L T = 45cm, and we x K i = 0.003S.cm -1 , K e = 0.003S.cm -1 and K T = 0.009S.cm -1 , and we choose the initial data as follows:

ρ(x, 0) = ρ min , U (x, 0) = 1 (ρ max -ρ min ) 2 .
(133) Moreover, in order to model electrical excitation within the myocardium (due to heart beat or external applied current), we apply a given external stimulus during a short time scale (0 < t -T d < T act , where T d is a time-delay), and we dene f is as (x = (x, y))

f is (x, t) = I app χ H (x)χ [0,Tact] (t -T d )χ prop (x, t -T d )Φ(x), (134)
where I app is the amplitude of the external applied stimulus with I app = 10 4 , and the functions χ H , χ [0,Tact] and χ prop are dened by

χ H (x) = 1 if x ∈ Ω H , 0 else (135) χ [0,Tact] (t) = 1 if t < T act , 0 else (136) χ prop (x, t) = 1 if x + y < 2t/T act , 0 else. ( 137
)
The last characteristic is associated to wave propagation as a diagonal which evolves from left bottom corner to right top corner. The function Φ, which corresponds to the shape of the electrical wave, is dened by

Φ(x) = 1 - 1 L H x - L H 2 2 - 1 L H y - L H 2 2 . ( 138
)
Figure 12 shows how the transmembrane potential ρ taken in the middle point of domain Ω H evolves according to the AP (after electrical wave passage), and we recover clearly the four phases. We present on Figure 13 the ionic behavior in this middle point with evolution of ionic variable U .

However, we can see these main phases with 3-dimensional representation of ρ and U at dierent moments. In a rst time, at t = 0.3s, the AP involves depolarization of the left bottom corner (Figure 14) of Ω H and ionic variable U begins to fall down (Figure 15). Then, at t = 0.5s (Figure 16), the electrical wave has covered more than the half of heart domain and we can see the dierent phases of AP: the left bottom corner is already on his recovery period (Phase IV), whereas its neighborhood is on the repolarization phase (Phase III) and we observe the plateau inherent to excited phase (Phase II) just behind the depolarized area (Phase I). In term of ionic variable, we also see the four phases with the activation of ionic channels in left bottom corner at t = 0.3s (Figure 15). At t = 0.5s (Figure 17), as for transmembrane potential, we can identify the four phases during the propagation of the electrical wave through domain Ω H .

Figures 18 and20 show that these four phases are visible through the extracellular potential ρ e at respectively t = 0.3s and t = 0.5s. Then, as Figures 24 to 27 show for ∆x = 1/30 to 1/100 in middle points of Ω H and Ω T , there is a few perturbation in the numerical solution, when a sudden change arises, due to chosen stopping criteria. This noise is tending to disappear totally with thinner discretization. Finally, we present on Figures 19 and21 the diusion of the electrical wave in the torso at t = 0.3s and t = 0.5s.

To complete the analysis, we study the relationship between the transmembrane potential ρ and ionic variable U , and their mutual inuence. For this, we introduce the following normalized terms

ρ N = ρ -ρ min |ρ max -ρ min | , U N = U -U min |U max -U min | . ( 139
)
We can see the correlation between transmembrane potential and ionic movements through the cell membrane in Figure 22. In particular, we clearly see that ionic variable begins to converge slowly to its equilibrium state at the beginning of the recovery period whereas transmembrane potential is already close to its equilibrium. Finally, in Figure 23 we perform the phase plan curve of ρ N and U N . After initialization at t = 0ms, the point A corresponds to equilibrium state (i.e. ρ = ρ min and U = 1/(ρ max -ρ min ) 2 ). We observe that ρ grows to ρ max in a shorter time scale before the ionic variable U starts to fall down (because ρ > ρ gate ). This depolarization phase (Phase I) leads us to point B. Then, between points B and C, we are in the excited phase (Phase II) which numerically corresponds to I ion (ρ, U ) = 0. After point C, ρ goes fast to ρ min while U varies slightly (Phase III). Finally, from point D to the end, ρ is near to its resting state and U goes back to its equilibrium state (Phase IV).

These analyzes led us to verify the eciency of our coupled modied LBM. Indeed, realistic behavior of the interactions between electrical potentials and ionic movements in heart-torso coupling (which is a critical numerical issue), has been demonstrated by the dierent simulations. The four phases of the AP are clearly reproduced. Our method remained stable and robust even if the solutions of elliptic equations contain some very local noises due to stopping criteria and sudden changes in dynamics.

Conclusion and commentary

An ecient and stable coupled LBM to solve a two-dimensional heart-torso coupled bidomain model is developed. From the Chapman-Enskog expansion analysis, the bidomain system which is a coupled of reaction-diusion, elliptic and ODE equations, can be correctly recovered by our modied LBM. This method is easy to implement. The preliminary results presented in this work show the capability of this method to capture the distribution of the electrical wave. It is clear that, due to the multi-scale nature of the system, the Cartesian grid used in our preliminary simulations is not very sucient to compute in a computationally ecient manner real life clinical situations (see e.g. [START_REF] Ellenbogen | Clinical Cardiac Pacing, Debrillation and Resynchronization Therapy[END_REF], [START_REF] Schulze | ECG Imaging of Ventricular Activity in Clinical Applications[END_REF] ) with complex geometry which is in general computationally expensive. Therefore, it is expected to solve the Lattice Boltzmann system on adapted Cartesian or triangular unstructured grids as e.g., in [START_REF] Valero-Lara | A Non-uniform Staggered Cartesian Grid Approach for Lattice-Boltzmann Method[END_REF][START_REF] Zhao | Adapted Unstructured LBM for Flow Simulation on Curved Surfaces[END_REF] and the references therein. Moreover, in order to overcome the limitations of the constraint CFL stability condition, we extend the method to implicit or semi-implicit time schemes, e.g., by using the θ-method (with θ ∈ [0, 1]) or Runge-Kutta methods, coupled with adaptive time stepping strategies, as e.g. in [START_REF] Huang | A Fully Implicit Method for Lattice Boltzmann Equations[END_REF] and the references therein. This coupled LBM method will be shown in a forthcoming paper for more general coupled models with realistic complex geometries. It would be interesting to use this developed method with observations coming from experimental data and a more complete description of the biophysical model of electrical cardiac activity. Moreover, since time delays in signal transmission are inevitable and a small delay can aect considerably the resulting electrical activity in heart and thus the cardiac disorders therapeutic treatment, it is then necessary to analyze the impact of delays on the dynamical behaviors of such a system by considering time-delays models with continuous and impulsive treatment (see e.g. [START_REF] Belmiloudi | Dynamical behavior of nonlinear impulsive abstract partial dierential equations on networks with multiple time-varying delays and mixed boundary conditions involving time-varying delays[END_REF]). In order to get even closer to a more realistic calculation, it is necessary to study, in the future, this method coupled with optimization technique and, optimal and robust control problems by using the approach developed in [START_REF] Belmiloudi | Mathematical modeling and optimal control problems in brain tumor-targeted drug delivery strategies[END_REF][START_REF] Belmiloudi | Stabilization, optimal and robust control. Theory and applications in biological and physical sciences[END_REF]. We have nondimensionalized ( 4)-( 6) by introducing the following dimensionless quantities (dimensionless variables, parameters and physical domains are denoted with " ∼ ") with ρ = ρ max -ρ min the characteristic action potential amplitude. The electrophysiological ionic state variable U is rescaled using

Ũ = ρ 2 U. (141)
Using these adimensional relations, we get nally the dimensionless system

Implementation and Algorithm

Our goal is to develop an ecient LBM for bidomain system. In this work, in order to validate the method, we just focus on stability and accuracy of the method for coarse discretization without any specic optimization. So we make our simulation in a beta version of our program, and the characteristic execution time can be signicantly reduce with using some Matlab functions or parallelizations. We want to control all the stages of build, and we need to represent every functions at each stage. That is why some stages are not optimized yet, even if we can do it easily (and so reduce the CPU time). The other main reason is that some stages are already implemented for complex boundary and complex domain for further works; the initialization stage or boundary treatment could be simplied for simpler boundaries. By the same way, we introduce some functions in order to have easy graphical representations. For all these reasons, we insist on the non optimal nature of our present code, and the characteristic CPU times we obtain are not revealing of the eciency of LBM. Moreover, the way we suggested for parallelization has not already been studied further by ourselves. Indeed, the two-stage build of the method guarantees that we can compute all local stages of the method with dividing the mesh between dierent threads. That includes local summations to recompose all the macroscopic values, the local denition of the equilibrium functions, second terms and link functions, and the local treatment of collision stage. Then, the streaming process can be resume as an index substitution between distribution function of adjacent nodes witch can also be separate for each velocity e i . This proposed parallelization will be done in a future work.

All the steps of this algorithm are represented in Figure 28. The step S1 is the initialization of the mesh (including vectors e i , vectors of boundaries index, or usual objects for graphical representations), the non-dimentionalization on time and space according to system (142), and denition of initial conditions. The step S2 is the treatment of the parabolic equation and of the ODE of system (142), and the both subroutines S3 and S4 correspond to the treatment of the elliptic equation of the system (142) in the heart and elliptic equation (143) in the thorax.

To estimate typical computionnal times, we have made our simulation with Matlab 2016b and with Intel core i5 processors. We present in Table 3 the dierent typical computational times for each part S1 to S4 of the algorithm in Ω × (0, T ) with T = 1 according to scale ∆x = 1/50, . . . , 1/100, ∆t = ∆x 2 , and convergence criteria parameters C ρe = C ρ T = 1 and N e = N T = 3. The execution time of both elliptic equations treatments are strongly related to the convergence criteria. With this simple implementation without optimization or parallelization, the execution time grows quickly according to the number of iterations and to the number of nodes. The rst scaling ∆x = 1/50 induces 2500 iterations to reach T = 1 and each domain contains 2601 nodes, and the last one ∆x = 1/100 induces 10 4 iterations for 10201 nodes for each domain. So we can conclude that execution time for each iteration is proportional to the number of nodes.

Figure 1 :

 1 Figure 1: Modeling of the membrane as resistor and capacitor coupled in parallel.

1 Figure 2 :

 12 Figure 2: Operator Hion and regularized operator Hr for rg = 1, rg = 0.5, rg = 0.1.

w 1 =

 1 w 2 = w 3 = w 4 = 1/9 and w 5 = w 6 = w 7 = w 8 = 1/36.

Figure 3 :

 3 Figure 3: Particle velocities for D2Q9 LBM Figure 4: Streaming process of a lattice node

κ

  200 (cm -1 ) transmembrane capacitance C m 10 -3 (F/cm 2 )Depolarization length

Figure 5 :

 5 Figure 5: Heart-torso domain.

Figure 6 :

 6 Figure 6: Potential ρ in ΩH at t = 1.

Figure 7 : 1 0Figure 8 :

 718 Figure 7: Ionic variable U in ΩH at t = 1

Figure 9 :

 9 Figure 8: Potential ρe in ΩH at t = 1.Figure 9: Potential ρT in ΩT at t = 1

  0.0937 7.972 × 10 -4 1/50 0.0334 2.039 × 10 -4 1/75 0.0141 7.971 × 10 -5 1/100 0.0074 4.483 × 10 -5 1/150 0.0032 1.990 × 10 -5 1/200 0.0017 1.121 × 10 -5 0.0541 3.117 × 10 -4 1/75 0.0277 1.471 × 10 -4 1/100 0.0187 7.314 × 10 -5 1/150 0.0097 3.417 × 10 -5 1/200 0.0053 2.237 × 10 -5

)Figure 11 :

 11 Figure 11: Error curves with respect to spatial step ∆x for U in ΩT at t = 0.5 and t = 1.

Figure 12 :

 12 Figure 12: Evolution of transmembrane potential ρ in middle of the heart (with four phases of AP).

Figure 13 :Figure 14 :Figure 15 :Figure 16 :

 13141516 Figure 13: Evolution of ionic variable U in middle of the heart.

Figure 17 :Figure 18 :

 1718 Figure 17: Ionic variable U in ΩH at t = 0.5s

Figure 19 :

 19 Figure 19: Potential ρT in ΩT at t = 0.3s

Figure 20 :

 20 Figure 20: Potential ρe in ΩH at t = 0.5s.

Figure 21 :

 21 Figure 21: Potential ρT in ΩT at t = 0.5s

Figure 22 :

 22 Figure 22: Evolution of ρN and UN in the middle of heart.

Figure 23 :Figure 24 :

 2324 Figure 23: (ρN , UN ) phase plan portrait.

Figure 25 :Figure 26 :

 2526 Figure 25: ρe for x = (7.5, 7.5) for ∆x = 1/30, 1/50, 1/75 and 1/100.

Figure 27 :

 27 Figure 27: U for x = (7.5, 7.5) for ∆x = 1/30, 1/50, 1/75 and 1/100.

  L∇, div = Ldiv, ∂/∂ t = T ∂/∂t,K ec = σ Kec , K ic = σ Kic , K T = σ KT ,where L is a characteristic length of physical domain Ω, T is a characteristic time and σ is a characteristic conductivity. The dimensional transmembrane and extracellular potentials ρ and ρ e , and thorax potential ρ T are rescaled using ρ = ρρ + ρ min , ρ e = ρρ e + ρ min , ρ T = ρρ T + ρ min(140)

  (.; ρ, Ũ ) -R h d iv( Kic ∇ρ) = R h d iv( Kic ∇ρ e ) + fis , on QH -R h d iv(( Kec + Kic ) ∇ρ e ) = R h d iv( Kic ∇ρ) + ( fes + fis ), on QH ∂ Ũ ∂ t = Hion (.; ρ, Ũ ), on QH(142)and -R h d iv( KT ∇ρ T ) = fT , on QT , (143) whereR h = σT L 2 κc m fis = T c m ρ f is , fes = T c m ρ f es , fT = T κc m ρ f T Ĩion (x, t; ρ, Ũ ) = T c m ρ I ion (x, t; ρ, U ),Hion (x, t; ρ, Ũ ) = T ρ 2 H ion (x, t; ρ, U ).

1 .

 1 If we take for exampleσ = L 2 κcm T , T = T and L = L H then R h = 1, QH = [0, 1] × [0, 1] × [0, 1] and Q = Ω × [0, 1], where Ω = [-L T L H , L T L H ] × [-L T L H , L T L H ].

Table 1 :

 1 Cell membrane parameters.

Table 3 :

 3 Typical computational time (in second) for each part of the algorithm in Ω × (0, 1).

	1/50	1.0598	58.9528	71.9957	217.8366
	1/75	7.8410 201.9650 243.5739	657.5601
	1/100 22.5423 764.5007 1042.5311 2567.0278

Acknowledgments: The autors would like to thanks all the referees for their relevant questions, suggestions and comments which lead to improve this paper.

Mesh definition, Non-dimentionalization, and initialization of initial conditions

Compute equilibrium fuction, second term and link function.

Collision and streaming steps.

Boundary treatment.

Compute ρ(t+∆t) and U(t+∆t).

Compute equilibrium function, second term and link function.

Collision and streaming steps.

Boundary treatment.

Compute ρe(t+∆t).

If convergence criteria is reached.