Kunkun, Boris, Dominik, Cyrile Gianluca

Edon Jean Marie

Maria

Meropi

Luka

Natalija

Przemek Thomas

Lorijn Sebastian Aga

Antonije, Pero Aleksandre

Keywords: électrophysiologie cardiaque, modèle bidomaine, modélisation multi-échelle, analyse asymptotique, homogénéisation, convergence à deux échelles, modélisation basée sur des images, jonctions communicantes, simulations numériques cardiac electrophysiology, bidomain model, multi-scale modelling, asymptotic analysis, homogenization, two-scale convergence, image-based modelling, gap junctions, numerical simulations

Montenegro et pour être toi-meme :) And to the more beautiful part of Carmen: Emanuela and Judit, my latin chatty girls, for all the intimate discussions that every girl needs, to Pauline and finally to Gwladys, pour être toujours lá pour moi depuis le premier jour, pour tout le soutien que tu m'as apporté et enfin pour m'avoir appris comment être the real untouchable.

I would like to say thanks to Manao team, especially to Xavier, Romain, and Nicolas who have helped me on the very beginning of my stay, and who understood the difficulty of the newcommer in Bordeaux. Further, I would like to thank to Geostat team and CQFD team for sharing the open space with Carmen and supporting a little bit of noise from time to time. :) I want to thank to all-boy-one-girl club known as Hiepacs team for making SIAM conference real fun, and for sharing beers and evenings afterwards in Bordeaux.

There have been many people who have passed through my life during my thesis and made my stay in Bordeaux more than just a work. Who were there to help me move 3 times, who helped me discover French way of life and many different cultures. First of all I have to thank Anna who was there from the first day and for staying my friend after all challenges I've been through.

Acknowledgement

I would like to thank to my adviser Yves Coudière for providing me with the opportunity to work on this PhD thesis, for teaching me about the cardiac electrophysiology and bidomain model and for giving me useful instructions during my struggle with the mathematical proofs and simulations. Also, for enabling me to participate in several international workshops, conferences and summer school. This gave me an exceptional overview of the current state of science and community in cardiology modelling. And of course for creating the great Carmen team (see below).

I would like to thank to my co-adviser Clair Poignard, firstly for accepting to join the project 5 months after it started, and then for teaching me about homogenisation and semigroup theory. I would like to thank him for dealing with my problems in full details and helping me to overcome all of the mathematical difficulties.

Then, I would like to give a very special thanks to Yves Bourgault, who worked with me intensively during my stay of three weeks in Ottawa University. This was a very inspiring, dynamical and fruitful work. I would like to thank him for making my stay in Canada very comfortable, and to his team for welcoming me as one of their own. I would like to thank to Thomas Desplantez for a very interesting collaboration, for giving his time to teach me all the details about gap junctions and for patience to reply my questions even if asked several times.

I would like to thank to Stephen Gilbert for providing me with data on the rat heart, and for being available for long and very useful discussions about the tissue structure and imaging.

I would also like to thank to Carmen members Jacques Henry, Mark Potse and Florian Caro for a very useful discussions about my work. And of course to Nejib Zemzemi, not only for useful advices for my work, but also for being a pleasant company during CINC conferences.

I would like to thank to my awesome officemates who made a life in open space fun and cheerful even during the most difficult of times: to Simon for being patient with my hundreds of questions per day on the very beginning of my thesis, to Mehdi for letting me use his hammock, for friendship but not for zombo.com, to Antoine, Julien, Cesare and Dragan, for sharing an open space, for organising parties and sharing drinks and dances, to Charlie for being the wise friend, for sharing with me a lot of Cameroon stories and for sharing the late evening office hours, to PEB for all nerdy discussions since our time together in Utah. And the last but not least to Marc, merci pour ton soutien pour la fin de ma thèse, pour m'avoir faite pratiquer le francais chaque jour, pour être venu au

Résumé

Dans cette thèse, nous avons abordé deux problèmes de modélisation mathématique pour la propagation des signaux électriques cardiaques : la propagation à l'échelle tissulaire en présence d'hétérogénéités et la propagation à l'échelle cellulaire avec des jonctions communicantes non linéaires.

Inclusions diffusives. Le modèle standard utilisé en électrocardiologie est le modèle bidomaine. Il est déduit par homogénéisation des propriétés microscopiques du tissu. Pour cela, on suppose que les myocytes électriquement actifs sont uniformément répartis dans le coeur. Bien que ce soit une hypothèse raisonnable pour des coeurs sains, ce n'est plus vrai dans certains cas pathologiques où des changements importants dans la structure tissulaire se produisent. C'est le cas, par exemple des maladies cardiaques ischémiques, rhumatismales et inflammatoires, de l'hypertrophie ou de l'infarctus. Ces hétérogénéités tissulaires sont souvent prises en compte à l'aide d'un ajustement ad hoc des paramètres du modèle. Le premier objectif de cette thèse consistait à généraliser les équations du modèle bidomaine au cas des pathologies cardiaques structurelles.

Nous avons supposé une alternance périodique d'éléments de tissus sains (modèle bidomaine) et modifiées (inclusions diffusives). La simulation numérique directe d'un tel modèle nécessite une discrétisation très fine, et entraîne un coût de calcul élevé. Pour éviter cela, nous avons construit un modèle homogénéisé à l'échelle macroscopique en utilisant une analyse à deux échelles. Nous avons retrouvé un modèle de type bidomaine avec des coefficients de conductivité modifiés, dits effectifs. En complément, nous avons effectué une vérification numérique de la convergence du modèle microscopique vers celui homogénéisé, dans une situation bidimensionnelle.

Dans la deuxième partie, nous avons quantifié les effets de différentes formes d'inclusions diffusives sur les coefficients de conductivité effectifs et leur anisotropie en 2D et 3D. De plus, nous avons effectué des simulations sur des domaines représentant des morceaux de tissu 2D avec ces coefficients de conductivité modifiés. Nous avons observé des changements de la vitesse de propagation et de la forme du front de l'onde de dépolarisation.

List of Figures

1.5

Membrane: the impermeable lipid bilayer is an insulator element that separate two different ionic solutions. However, it is permeated by the presence of ionic channels, producing a potential difference. From [START_REF] Rudy | [END_REF] . . . . . . . . .

1.6

Action potential of a ventricular muscle cell. From [START_REF] Collin | Asymptotic analysis in cardiac electrophysiology: applications in modeling and in data assimilation[END_REF]. . . . . . . . . . . .

1.7

The equivalent electrical circuit for a membrane with three types of ionic channels and a leakage current. Created by Nrets, CC BY-SA 3.0, https: //commons.wikimedia.org/w/index.php?curid=17095070. . . . . . . . . .

1.8

Microscopic structure used to derive the macroscopic, bidomain model. . . .

1.9

Reconstruction of the fiber orientation in the rabbit heart using DT-MRI. (A) Lon-axis cross section showing fiber helix angle. (B) Transmural fiber helix angles extracted from a sector of the left ventricular wall. Taken from [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF].

1.10 Details of the gap junction: each channel is made of two hemichannels or connexons. Each connexon is running transverse to the opposing membrane. Connexons of the adjacent cells align themselves to form the structural continuity. Edited image taken from [START_REF] Makowski | Gap junction structures: Analysis of the x-ray diffraction data[END_REF]. . . . . . . . . . . . . . . . . . . . . . . 1.11 Configurations of connexons and gap junction channels for two different connexins (Cx). Fourteen different connexons can be formed when two connexins are co-expressed in a cell: 2 homomeric connexons, consisting of one type of connexin, and 12 heteromeric connexons consisting of two types of connexins. Gap junction channels formed after docking of two connexons can be homotypic when the connexons are identical (composed of one or two types of connexins) or heterotypic when the connexons are different (composed of one or two types of connexins). Taken from [START_REF] Desplantez | Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties[END_REF]. . . . . . . . xvii 1.12 Left: Confoncal microscopic projection image of an optical section series through a single isolated ventriculat myocyte labeled for Cx43, showing the positions of clusters of transversely oriented GJCs (arrows). Scale bar = 20 µm. Taken from [START_REF] Peters | Reduced content of connexin43 gap junctions in ventricular myocardium from hypertrophied and ischemic human hearts[END_REF]. Right: Lattice of regular hexagons, representing the transverse-section of GJ used in the model in [START_REF] Wilders | Limitations of the dual voltage clamp method in assaying conductance and kinetics of gap junction channels[END_REF]. The diameter of a gap junction is d, the diameter of a gap junction channel is d c ≈ 2 nm, and the lattice constant is d hex ≈ 9 nm. . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Left: DIC image (Differential Interference Contrast) of a pair of cells imaged during dual whole cell patch clamp experiment. Taken from [START_REF] Del Corsso | Transfection of mammalian cells with connexins and measurement of voltage sensitivity of their gap junctions[END_REF]. Right: The equivalent electrical circuit of the dual whole cell patch clamp experiment. .

2.1

On the left: the idealised full 2D domain, Ω. On the right: the periodic cell, Y .

2.2

Plots for v ε , for different size of ε, and for the homogenised v. We can observe the convergence of the solutions. . . . . . . . . . . . . . . . . . . . . . . . .

2.3

Convergence study for v ε and h ε . Observed convergence rate for the given time: 1.39 for V , and 0.63 for h. . . . . . . . . . . . . . . . . . . . . . . . . .

3.1

Unit cell geometries, symmetric w.r.t. x-and y-axis, used to test the effects of the diffusive inclusions on the conductivity tensors in 2D. Case I: the shape of the inclusion is a circle. The volume of the inclusion can go from 0% to 76% of the unit cell. Case II and III: the shape of the inclusion is an ellipse with r x : r y = 2 : 1 and 1 : 2, resp. The volume of the inclusion can take up to 39% of the unit cell. Case IV: the shape of the inclusion is a superellipse with e = 0.2 and r x : r y = 1 : 1 . The volume of the inclusion can take up to 98% of the unit cell. Case V and VI: the shape of the inclusion is a superellipse with e = 0.2, and r x : r y = 2 : 1 and 1 : 2, resp . The volume of the inclusion can take up to 49% of the unit cell. . . . . . . . . . . . . . . . . . . . . . . .

3.2

Unit cell with the nonsymmetric inclusion. The volume of the inclusion can take up to 63% of the unit cell. . . . . . . . . . . . . . . . . . . . . . . . . .

3.38

Initial condition for the simulation where we try to obtain the spiral. The domain size is [-30, 30] × [-30 , 30]. On the left is the scaled transmembrane voltage v 0 , which is zero everywhere except in the small region -30 < x < -17 and -3 < y < 3. On the right is the gating variable h 0 that takes a constant value 0 under x-axis, and the value h 0 = 0.5 above x-axis. . . . . .

3.39

Simulation of the standard bidomain model. The tissue is assumed to be homogeneous, with constant conductivities. We observe the beginning of the spiral, but it quickly disappears. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms, ... , 480ms.

SCARS:

The heterogeneous structure of the tissue. We plot the volume fraction |Y B | that ranges from 0.0 for the "dead" center of the scar, to 1.0 for the healthy tissue. The rings where the value ranges 0 < |Y B | < 1, we call the border zone. On the left (SCAR L) the center of the scar has radius r 1 = 16, in the middle (SCAR M) r 1 = 9, on the right (SCAR S) r 1 = 4. For all tree cases the whole scar has the radius r 2 = 22. . . . . . . . . . . . . . . . . . . 3.42 SCAR L. Case IV. Simulation of the modified bidomain model, with the inclusions assumed to be of the type IV. The tissue is assumed to be heterogeneous as in SCAR L. The scar slows down the propagation of AP in the center of the tissue, and forces the wavefront to go around the scar. This provides enough time for the tissue to go through the recovery period, and to be reactivated. Hence, we observe the spiral that is guided by the shape of the scar. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms, ... , 1120ms. . . . . . . . . . . . . . . . . . . . .

SCAR M.

Case IV. Simulation of the modified bidomain model, with the inclusions assumed to be of the type IV. The tissue is assumed to be heterogeneous as in SCAR M. Like before, the scar is big enough to stop the propagation of AP in the center of the tissue, and forces the wavefront to go around the scar. Hence, we observe again the spiral that is guided by the center of the scar. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 60ms, 120ms, 180ms, ... , 1680ms. . . . . . . . . . . . .

3.44 SCAR S. Case IV. Simulation of the modified bidomain model, with the inclusions assumed to be of the type IV. The tissue is assumed to be heterogeneous as in SCAR S. In this case the center of the scar is small, and not enough to stop or slow down the propagation of AP trough the center. Although, we do observe a change in the shape of the wavefront, we lose the spiral. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms, ... , 480ms. . . . . . . . . . . . . . . . . . . . . .

3

.45 SCAR M. Case VI. Simulation of the modified bidomain model, with the inclusions assumed to be of the type VI. The tissue is assumed to be heterogeneous as in SCAR M. There is a change of the direction of the propagation of AP in the scar, which alternates the shape of the wavefront significantly. This kind of activation does not provide the spiral in the tissue in contrast to the same setting with the inclusions of the type IV.

4.9

The bidirectional vector field of the fiber orientations. Data obtained from Giblert et al. 2012 [START_REF] Gilbert | Visualization and quantification of whole rat heart laminar structure using high-spatial resolution contrast-enhanced MRI[END_REF], the original data on the full heart (left) and the projection of the fibers on the segmentation (right). Visualised with SCIRun.

4.10

The bidirectional vector field of the fiber orientations projected on the mesh nodes, visualised with SCIRun. This field is much finer than the one in the Figure 4.9, due to the much higher number of nodes in the mesh, compared to the number of voxels in DT-MR image.

5.5

Transjunctional current for Cx43 simulated with the parameters obtained from the experiments and the data fitting + corrected G j,max . . . . . . . . .

5.6

Geometric settings for 2D mathematical model. Ω i , Ω e -intra and extracellular domains, Γ ion -"ionic" boundary, Γ j -"junctional" boundary. . . . . . . . . . . . .

5.7

Preliminary numerical results on 4 myocytes. . . . . . . . . . . . . . . . . .

5.8

Preliminary numerical results for the microscopic model without gap junctions. We do not observe propagation of AP from one cell to the other. . . . . . . .

List of Tables

1.1

Parameters of the bidomain model for a human heart. Three different dimensional scalings are proposed. Adim1: based on the transmembrane potential upstroke. Adim2: based on the action potential duration and width. Adim3: based on convenient units. The table is taken from [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF]. . . . . . . . . . . .

3.1

The bidomain conductivities for ventricular tissue, values in [S m -1 ]. The first 5 columns of the table taken from [START_REF] Sundnes | Computing the electrical activity in the heart[END_REF], and the last column from [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF].

3.2

The test cases for geometries of inclusions in 2D. . . . . . . . . . . . . . . .

3.3

The test cases for the geometries of the inclusions in 3D. . . . . . . . . . . .

Heart

This study started since more than two millennia back. In the time of Aristotle (350 B.C.) the heart was considered to be the most important organ in the body, while other organs like the brain and lungs, were thought to exist to cool the blood. Until this day the heart keeps the position of one of the most important and the most studied organs in the human body.

Physiology and function

The heart is located in the lungs, it has a size of the fist and weights 200-300 g. The heart is composed of two separate pumps, see Figure 1.1. The right heart pumps the

Electrophysiology

The most important property of cardiac cells in this sense is that they are excitable, i.e. they have the ability to respond actively to electrical stimulus. In fact, the electrical activity is a consequence of the chemical processes and exchange of electrically charged ions between the interior of the cells and it's surrounding.

In the resting state the cell maintains internal ionic concentrations different from those outside of the cell. Hence, there is the difference of electrical potential across the cell membrane, which is called the transmembrane voltage. When an electrical stimulus is applied to a cell there is a change in the transmembrane voltage. If the magnitude of this change is above a certain threshold value, the conductive properties of the cell membrane change and there is a rapid flux of ions through the cell membrane. This in turn causes a rapid change in the transmembrane voltage, from the -85 mV in the resting state to about +20 mV. This process is called depolarisation, and happens in about 0.2 s in the human heart, while in the other animals it can happen in less than 10 ms. After depolarisation, the cardiac cell goes through a "long" plateau phase and then it returns to the resting state, i.e. it repolarises. The complete cycle of depolarisation, plateau phase and repolarisation is called the action potential (AP).

The cardiac cells are connected by gap junctions, intercellular pathways made of proteins, that are permeable to ions and small molecules. They provide a direct electrical connection between the neighbouring cells. When one cell is depolarised this connection will affect the potential in the neighbouring cells, and if the change is above a threshold value it will trigger the neighbouring cell to go through it's own AP. In other words, gap junctions enable the action potential propagation though the heart.

A propagation of the AP in the whole heart generates the rhythmical heart beat called the sinus rhythm. The main elements of the cardiac electrical conduction system are shown on the Figure 1.2. The electrical signal in the heart starts in the sinoatrial node (SA node) which is located above the right atrium. The cells of the SA node are self-oscillatory cells, called pacemaker cells. They spontaneously produce the action potential and do not need external stimuli.

The activation of a small region in the atria by pacemaker cells, results in a wavefront of depolarisation, which activates the complete atria and causes them to contract. The atria are separated from the ventricles by a non-conducting tissue called septum. Hence, the only way to pass the electrical signal to the ventricles is through the atrioventricular node (AV node), found at the bottom of the right atrium. Due to the small conduction velocity in the AV node there is a small delay between the activation of atria and ventricles, which enables a better filling of the ventricles and improves the pumping function of the heart.

The signal passes from the AV node to the bundle of His, that branches out into the left and right bundle. In the ventricles the signal is carried by a specialized tree-like structure, called the Purkinje fibers, which provide a rapid conduction of the electrical signal. The ventricular muscle cells are stimulated at the ends of Purkinje network, which are located near the inner surface of the ventricles, called endocardium. This results again in a wavefront of depolarisation that activates both ventricles and makes them contract. To find out more details on the heart physiology one can look into [START_REF] Hall | Textbook of medical physiology[END_REF].

The electrical activity is an essential function of the heart, and many heart diseases are closely linked to disturbances in the electrical activity: heart failure, cardiomyopathy, heart murmurs, etc. Even though it is a well studied problem and there is already a vast literature on the topic, there are still many mechanisms that are not fully understood. For example, one of the biggest mysteries is the ventricular fibrillation. It is the state in which cells contract in a seemingly random manner, so the heart is unable to pump the blood.

The problem in understanding all the mechanisms in the heart come from the complexity of combining all the processes on several scales. The organ-level electrical activity is the result of the billions of small scale processes occurring in the cells, and understanding of how all these processes interact is a very difficult task. Mathematical modelling and computer simulation are the new techniques that help us improve our understanding of the large scale processes, but also they save us the cost and time of conducting real-life experiments, see the flowchart in the Figure 1.8. Further in this chapter we give a short overview of the mathematical models commonly used to describe the electrical activity of the heart, from a cell scale to the organ and body scale, and point out some limitations of these models that will be addressed in this thesis.

Pathology statistics

Before we dive into details about the mathematical modelling of the cardiac electrophysiology, let us give a motivation to study the heart into details.

Cardiovascular diseases (CVD) are a major health problem in Europe and in the industrialised world. According to the 2012 European Cardiovascular Disease Statistics [START_REF] Nichols | European cardiovascular disease statistics[END_REF], published on http://www.escardio.org/, each year CVD causes over 4 million deaths in Europe and 1.9 million deaths in the European Union (EU). That makes 47% of all deaths in Europe and 40% in the EU. CVD are the main cause of death in women in all countries of Europe and the main cause of death in men in all but 6 countries. See Figure 1.4. Overall CVD are estimated to cost the EU economy almost C196 billion a year. Additionally, CVD are the main cause of death before the age of 75, and the second cause of death (after cancer) before the age of 65. While the medical practitioners are looking into the causes and correlations between CVD and diet, physical activity and a lifestyle, we are on the quest to understand the underlying biophysical mechanisms in cardiac physiology and pathology. An improved understanding of how the heart works might lead to new techniques for the diagnosis and treatment of heart problems.

Cell scale modelling

The cell is enclosed by a membrane whose thickness is about 7.5 -10.0 nm. The membrane consists of two layers of phosphoglycerides molecules. These molecules have hydrophilic heads and hydrophobic tails, which means that the molecule heads interact with the outside and the inside of the cell, and the tails are pointing inside of the membrane, see Figure 1.5.

All over the membrane, embedded in the lipid bilayer we can find the ionic channels. They are macromolecular pores through which sodium, potassium, chloride flow inside and outside of the cell. The flow of these ions forms the basis for the bioelectric phenomena, [START_REF] Malmivuo | Web-version of the book: Bioelectromagnetism -Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF]. The ionic channels are made of large proteins, that are usually very specific and permeable only to a one kind of substance. In mathematical models, it is necessary to include the flow of all the ions through these channels, hence many variables need to be used to characterize the dynamics of the channels as they open and close.

The story of the ionic models and bioelectrical modelling can not pass without mentioning a great breakthrough of Alan L. Hodgkin and Andrew F. Huxley in the 1950s. They have worked on the giant axon of the squid (a part of the nerve cell) and created a Fig. 1.5: Membrane: the impermeable lipid bilayer is an insulator element that separate two different ionic solutions. However, it is permeated by the presence of ionic channels, producing a potential difference. From [START_REF] Rudy | [END_REF] model for the signal propagation that was based on the detailed models for several ionic currents in the cell.

Cell membrane

Electrophysiologists model the effects of ionic concentration differences, ion channels, and membrane capacitance in terms of an equivalent electrical circuit. Ludvig Hermann (1905) was the first to suggest that under subthreshold conditions the cell membrane can be described by a uniformly distributed leakage resistance and parallel capacitance. The differential equation for this simple circuits reads,

C m dv dt = -I ion , (1.1)
where C m is the capacitance of the membrane per unit area, I ion is the ionic current per unit area, and v is the transmembrane voltage. The capacitance is determined by the properties of the lipid bilayer, and is taken to be fixed. In a simplified version I ion = Gv, where G is a fixed conductance of ionic channels. The solution of this equation for v is the exponential decay in time. Such a model is very simple to understand and solve but it gives only the starting point for real-life modelling. In real life the conductance of most of the ionic channels are not constant in time and are environment dependent. Changes in the conductance can happen in less than 1 ms, which depend primarily on the voltage. Hence, they are called time and voltage dependent channels. However, the differential equation used to model a membrane patch is commonly a modified version of a resistor-capacitor circuit equation.

A typical action potential for a ventricular muscle cell is given in the Figure 1.6. The value of the transmembrane potential in the resting state is around -85 mV in the human heart, while it can range from -70 mV to -100 mV in other species, [START_REF] Sundnes | Computing the electrical activity in the heart[END_REF]. Such a complicated AP is a consequence of many different kinds of ionic channels in the cell. In cardiac muscle, the action potential is caused by opening of two main types of channels: (1) the fast sodium channels and (2) another entirely different population of slow calcium channels, which are also called calcium-sodium channels. This second population of channels differs from the fast sodium channels in that they are slower to open and, even more important, remain open for several tenths of a second. During this time, a large quantity of both calcium and sodium ions flow through these channels to the interior of the cardiac muscle fiber, and this maintains a prolonged period of depolarization, causing the plateau in the action potential. The calcium ions that enter during this plateau phase activate the muscle contractile process, [START_REF] Hall | Textbook of medical physiology[END_REF].

Commonly, each ionic channel is described separately, with the equivalence circuit given in the Figure 1.7. The differential equation for this circuits is,

C m dv dt = -I N a -I K -I Cl -I Leak , (1.2)
Each of the ionic channels has different time and voltage dependent behaviour as will be explained in the following text.

Ion flux

The ionic channels can be passive, when the ionic concentration follows the electrochemical gradient. This is the case, for example, of the fast N a channel. Other ionic channels behave as pumps or exchangers, which are important for maintaining the correct ionic concentration in the cell. These are able to transport ions in the direction opposite to the electrochemical gradient. This process is accomplished either by using the concentration gradient of different ion or by consuming chemically stored energy in the form of ATP. This is the case, for example, of N a -K exchanger, [START_REF] Sundnes | Computing the electrical activity in the heart[END_REF].

It is not easy to assess the exact behaviour of the channel in the course of time. The first step in deriving a model for the ionic channel is observing the properties of the membrane in the equilibrium state. Considering only one type of ion one can derive the equilibrium equation, where the diffusive flux and the electrically driven flux are equal.

I ion = I D + I E = 0. (1.3)
The diffusive flux is given with the Fick's law:

I D = -D∇c, (1.4)
where D is the diffusion coefficient and c is the concentration of the ion. The electrically driven flux is given by Planck's equation:

I D = -m z |z| c∇φ, (1.5)
where z is the charge, φ the scalar electrical potential, and m is a mobility which is related to the diffusion coefficient as

m = D |z|F RT , (1.6)
where R is the ideal gas constant, T is the temperature, and F is Faraday's constant. Finally, the total flux of the ion gives the Nernst-Planck equation for one ion,

I ion = -D ∇c + zF RT c∇φ . (1.7)
In the equilibrium state of the cell, i.e. when the total ionic current is zero, there will be a potential difference across the membrane, which we refer to as a transmembrane voltage, v = φ i -φ e , where φ i and φ e are the intracellular and extracellular potentials. When we solve the equation (1.7) for the equilibrium state we find that the transmembrane potential that gives zero flux is given as

v eq = RT zF log c e c i , (1.8)
where c i and c e are the intracellular and extracellular ion concentrations. The voltage v eq is called Nernst equilibrium potential.

In the non-equilibrium the derivation of the actual flux is more complicated. We know that in any case, the flux must satisfy the condition to be zero when the transmembrane potential is v = v eq . The simplest model for the flux of a given ion following the electrochemical gradient is the linear model,

I ion = G(v -v eq ), (1.9)
where G is the electrical conductance for this ion. Depending on the ion of interest, G might be a constant, or a function of time, membrane potential or ionic concentration. There are other models for the ionic flux that satisfy Nernst equilibrium potential. The most commonly known is the Goldman-Hodgkin-Katz formulation. For further reading on these models one can look at [START_REF] Malmivuo | Web-version of the book: Bioelectromagnetism -Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF][START_REF] Sundnes | Computing the electrical activity in the heart[END_REF]. The next step is to understand what happens if the membrane is permeable to more than one ion. We can find the Nernst equilibrium potentials for each ion separately from (1.8),

v 1 eq = RT zF log c 1 e c 1 i , v 2 eq = RT zF log c 2 e c 2 i , • • • (1.10)
In general the equilibrium transmembrane potentials are not the same, hence it is not possible to achieve the equilibrium for all ions simultaneously. What we can expect to have is that the total current is zero,

I ion = I 1 + I 2 + • • • , (1.11)
and find the equilibrium transmembrane potential in this case. For example, if we deal with only two types of ions and use the linear model for each ionic flux, we have

G 1 (v -v 1 eq ) + G 2 (v -v 2 eq ) = 0. (1.12)
From here easily follows,

v eq = G 1 v 1 eq + G 2 v 2 eq G 1 + G 2 .
(1.13)

Gating

Ionic channels are described as series of gates, that open or close in response to an electrical stimulus. The mathematical description of the ionic channels considers the channels to be composed of several sub-units. Every sub-unit may be either open or closed, and ions can pass only when all units are open. If we consider a channel that contains only one sub-unit, the change between open and closed state may be written as

C α β O, ( 1.14)
where C stands for closed, O for open, and α, β stand for the rates of opening and closing, respectively. These rates are normally voltage dependent. Let say that the probability for the channel to be in the open state is g, then the dynamics of opening and closing is described as

dg dt = α(v)(1 -g) + β(v)g. (1.15)
The solution to this equation in general is,

g(t) = g 0 + t 0 α(v)e s 0 (α(v)+β(v))
dr ds e -t 0 (α(v)+β(v))ds .

(1.16)

If we assume that α and β are constant the equation can be rewritten as

dg dt = g ∞ -g τ g , ( 1.17)
where

g ∞ = α α + β
, and

τ g = 1 α + β . (1.18)
Then the solution is

g(t) = g ∞ + (g 0 -g ∞ )e -t/τg , (1.19)
where g 0 is the initial value for g, and g ∞ is the steady state value that g approaches as t increases. The rate at which g approaches g ∞ is determined by the time constant τ g . In reality α and β are voltage dependent and one should use the solution (1.16). On the other hand, if τ g is small enough, one can assume that g reaches the steady state value almost immediately, and g ∞ can be used as a reasonable approximation for g(t).

When the channel consists of more sub-units the probability that the whole channel is open depends on the probability of each sub-unit to be open. Assuming that the sub-units are an independent gates, i.e. they open and close independently one from another, then the probability that the whole channel is open is the product of each subunit being open. So if the channel consists of n identical sub-units, then the probability that the whole channel is open is given as

O = g n .
(1.20)

If the channel consists of m units of type g, and n units of type h, then the probability will be

O = g m h n , (1.21)
where the dynamics of the channels g and h is given in the form of the equation (1.15), with respective rates α g , β g and α h , β h . Finally, the total current through the ionic channels of the same kind can be described as a product of the maximum current possible, when all channels are open, and the proportion of all open channels. If we use the linear model as in (1.9), the ionic current is given as

I = G max O(v -v eq ), (1.22)
where G max is the maximum conductance if all channels are open and v e q is the equilibrium potential for the specific ion.

Ionic models

The ionic models are used to describe the dynamics in excitation of the cardiac cells, i.e. the cardiac action potential (AP). There are numerous ionic models developed to match different kinds of cardiac cells depending on the species, age, or even the position of the cells in the heart. There is a large database of different ionic models published to date on the website www.cellml.org. Ionic models are commonly grouped in three cathegories:

• the physiological cell ionic models (e.g. Luo-Rudy, Beeler-Reuter, Hodgkin-Huxley)

• the reduced ionic models (e.g. the 3-variables Fenton-Karma model is a reduction of the Luo-Rudy 1 model),

• the phenomenological models (e.g. FitzHugh-Nagumo, Aliev-Panfilov, Mitchell-Schaeffer).

The Hodgkin-Huxley model

The Hodgkin-Huxley model is historically the most important breakthrough in modelling the activity of excitable cells. They worked on the giant axon of the squid (not the giant squid), [START_REF] Keener | Mathematical Physiology I: Cellular Physiology, II: Systems Physiology[END_REF]. Their model was originally used to explain the AP in the long giant axon of a squid nerve cell, but the idea was since used to describe the AP of many other excitable cells. For the historical importance we give the explanation of their model here.

As we said earlier, the cell membrane can be modelled as a capacitor in parallel with several ionic currents. In the squid's giant axon the principal ionic currents are the N a + current and the K + current. In the Hodgkin-Huxley model all other currents are combined together into one current called the leakage current.

C m dv dt = -I N a -I K -I L + I app , (1.23)
where C m is the membrane conductance, and I app is the initial current applied in order to trigger the action potential. The instantaneous I -V curves for open N a + and K + channels are approximately linear, so the model reads,

C m dv dt = -g N a (v -v N a ) -g K (v -v K ) -g L (v -v L ) + I app , (1.24)
where g N a , g K , g L are the respective conductances. These dependences are voltage dependent and the main challenge for Hodgkin and Huxley and their collaborators at Cambridge was to estimate the conductances of each channel separately. Thanks to the voltage clamp technique they were able to keep the transmembrane voltage on a constant, known level and observe the ionic currents, i.e. the conductance in time. Furthermore, they were able to separate the total ionic current into its constituent parts. They have made an important observation that when the voltages are fixed, the conductances are time dependent. The rest was the mathematical modelling by fitting the experimental data, without the physiological meaning. They came up with the following expressions for the conductivities that gave acceptable agreement with experimental data.

g N a = G N a m 3 h, (1.25) g K = G K n 4 , (1.26) g L = G L , (1.27)
where G N a , G K , G L are constants that represent the maximum conductance for the corresponding channels, and m, h, n are voltage dependent variables that obey simple differential equations, such as (1.15). We call them the gating variables.

To sum up, the Hodgkin-Huxley model is given with the differential equation,

C m dv dt = -G N a m 3 h(v -v N a ) -G K n 4 (v -v K ) -G L (v -v L ) + I app , (1.28)
with the gating variables given by equations of the form,

dg dt = α g (v)(1 -g) -β g (v)g, (1.29)
for g = m, h, n, and α g , β g area the rate functions that depend on v. The rate functions are determined from the experimental data. For more details about this model one may look into [START_REF] Keener | Mathematical Physiology I: Cellular Physiology, II: Systems Physiology[END_REF][START_REF] Plonsey | Current flow patterns in two-dimensional anisotropic bisyncytia with normal and extreme conductivities[END_REF][START_REF] Sundnes | Computing the electrical activity in the heart[END_REF].

The FitzHugh-Nagumo model

The FitzHugh-Nagumo (FHN) model was suggested about 10 years after the Hodgkin-Huxley model. The idea is to extract the essential behaviour of the Hodgkin-Huxley phase plane and represent it in a simplified form. This was achieved elegantly, with only two variables: the fast one (v) and the slow one (w), also called the recovery variable. The FHN model formulation is given by

dv dt = c 1 v(v -a)(1 -v) -c 2 w + I app , (1.30) dw dt = b(v -c 3 w). (1.31)
Here a, b, c 1 , c 2 , c 3 are given parameters that can be adjusted to simulate different cell types. In the original model these values were set to approximate the Hodgkin-Huxley model for the AP of the giant axon. This model is able to reproduce the most important characteristics of the AP. Such simplified models are useful in simulations as they have a small number of variables and are easier to solve, but they are limited when one needs to study how the changes in physiology on a cellular level affects the function of tissues and complete organs.

The Noble model

The first model to describe the AP of the cardiac cell was proposed by D. Noble in 1962, and it was a model for the Purkinje cells. It was based on the Hodgkin-Huxley model with the parameters changed to fit the AP of the Purkinje cells. The main difference between the AP of the squid giant axon and the those of Purkinje cells is the duration of the plateau phase, that is 300-400 ms in cardiac cells compared to 3 ms in axon. Still, the model describes two ionic currents: potassium and sodium, with leakage current assumed to be zero. Additionally, Noble assumed two separate potassium currents I K1 and I K2 .

The model is given as,

C m dv dt = -g N a (v -v N a ) -(g K1 + g K2 )(v -v K ) + I app , (1.32)
where,

g N a = G N a m 3 h + G i , (1.33
)

g K1 = G K1 n 4 , (1.34
)

g K2 = g K2 (v). (1.35)
As we can see the expressions for the N a current and K1 current are similar to those of the Hodgkin-Huxley model, and they have time dependent conductances, while g K2 is time independent. The G i conductance, represents a leakage current of sodium that serves to prolong the AP. At the time when the model was first written there were no detailed data on the ionic currents in the cardiac cells. However, the model was a successful attempt to reproduce the AP of Purkinje cells. It was finally adjusted in 1975 to include more ionic currents and to give more accuracy on the underlying physiology.

The Beeler-Reuter model

Another important model for simulating AP in cardiac ventricular cells came in 1975, and it was proposed by Beeler and Reuter. The important addition in this model is the role of calcium, which is essential for the contraction of the muscle cells.

This model describes four ionic currents: the fast inward sodium current, a slow inward calcium current, and two outward potassium currents. These four currents are controlled with six gating variables. Additionally, the model describes the intracellular calcium concentration. The model reads,

C m dv dt = -I N a (m, h, j, v) -I K (v) -I x (x, v) -I s (f, d, v) + I app , (1.36) dc dt = 0.07(1 -c) -I s , (1.37)
where c represents the calcium concentration. The gating variables: m, h, j, x, f, d satisfy differential equations

dg dt = α g (v)(1 -g) -β g (v)g, (1.38)
where g = m, h, j, x, d, f , and as before α g , β g represent the respective rate functions.

The Beeler-Reuter model is still widely used in the computations. The reason for this is the fact that it has only eight state variables, while the more recent models are more and more complex and stiffer, including more ionic channels and many state variables, and are computationally much more challenging. This is an important factor in simulations of the tissue patches or the whole heart, which are computationally very demanding.

The Mitchell-Schaeffer model

The model that we mostly use in our current work is the Mitchell-Schaeffer (MS) model proposed in 2003, [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF]. It is a model for electrical activity of cardiac membrane which incorporates only an inward and an outward current. It is obtained as a special case of the Fenton-Karma model (1994), with a proper rescaling.

The voltage is dimensionless and rescaled so it ranges between zero and one, and is governed by the ODE

dv dt = I in (v, h) + I out (v) + I app , (1.39)
where the three currents are given as follows: I in is a combination of all currents which rise the voltage across the membrane (primarily sodium and calcium) and is given as.

I in (v, h) = 1 τ in hv 2 (1 -v), (1.40)
where h is the gating variable. This current mimics the behaviour of the fast gates, and is specified by the time constant τ in . The outward current is a combination of all currents that decrease the transmembrane voltage (mainly potassium), and is given as

I out (v) = - v τ out , ( 1.41)
with the time constant τ out . So the full model can be written as,

dv dt = 1 τ in hv 2 (1 -v) - v τ out + I app , (1.42) dh dt = g(v, h). (1.43)
The dimensionless gating variable h, varies between 0 (closed) and 1 (fully open), and satisfies the following ODE, This model is useful for its simplicity, it is comparable to the FitzHugh-Nagumo model, which makes it useful for numerical simulations, especially in two or three spatial dimensions where numerical efficiency is important. Several modified versions of the MS model may be found in the literature. We will make use of the one proposed by Boulakia et al. 2008, [13]. The model is a regularised version of the MS model, and it reads

dh dt =    1-h τopen if v ≤ v gate , -h τ close if v > v gate , ( 1
I ion = 1 τ in hv 2 (1 -v) - v τ out , (1.45) g(v, h) = 1 τ close + τ close -τ open τ close τ open h ∞ (v) (h -h ∞ (v)), (1.46)
where,

h ∞ (v) = 1 2 1 -tanh v -v gate µ gate , ( 1.47)
and µ gate > 0. This regularised version of the MS model comes in handy in the mathematical manipulations, as we get rid of the jump in the function g(v, h).

Tissue scale modelling

The cardiac tissue consists of the network of cardiac cells, each surrounded with the extracellular matrix and connected to the adjacent cells through the gap junctions. If one applies the electrical stimulus in one cell, it will go through the AP, and the electrical signal will get propagated to the neighbouring cells. This process is complicated to model even if we deal with a few cells. Knowing that the cardiac tissue contains millions of cardiac cells, such discrete modelling would be mathematically very challenging and numerically extremely expensive.

Extensive studies on the electrical properties of the biological tissues and especially on the cardiac tissue, were conducted in the 70's and 80's, [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF][START_REF] Geselowitz | A bidomain model for anisotropic cardiac muscle[END_REF][START_REF] Walter | Simulation studies of the electrocardiogram. I. The normal heart[END_REF][START_REF] Muler | Electrical properties of anisotropic neuromuscular syncytia. I. Distribution of the electrotonic potentia[END_REF][START_REF] Muler | Electrical properties of anisotropic neuromuscular syncytia. II. Distribution of a flat front of excitation[END_REF][START_REF] Muler | Electrical properties of anisotropic neuromuscular syncytia. III. Steady state of the front of excitation[END_REF][START_REF] Peskoff | Electric potential in cylindrical syncytia and muscle fibers[END_REF][START_REF] Peskoff | Electric potential in three-dimensional electrically syncytial tissues[END_REF][START_REF] Tung | A bi-domain model for describing ischemic myocardial dc potentials[END_REF]. What has been understood from the experiments is that the intracellular space behaves as if interconnected, but with higher conductivity along the fiber than transverse to it. The same holds for the extracellular, or interstitial space. Schmitt (1969) proposed the macroscopic model of the two domains that are considered to be continuous and interpenetrating. As we will explain later, the simple Ohm's law can be used to describe the electrical current in these two regions and the transmembrane current flowing from one to another region, to obtain what is commonly known as the bidomain model. Neu andKrassowska (1992, 1994) [START_REF] Neu | Homogenization of syncytial tissues[END_REF][START_REF] Krassowska | Effective boundary conditions for syncytial tissues[END_REF], suggested the first mathematical approach in understanding the derivation of bidomain model from the microscopic description of the tissue. The derivation is based on an idealized representation of the syncytium as a network of interconnected cells arranged periodically in space. From this microscopic model, a homogenization process derives a macroscopic volume-averaged model for intracellular, extracellular and transmembrane voltage and its boundary conditions. This model corresponds to the already used bidomain model, with deeper understanding of the parameters and their connection to the cell scale processes.

The bidomain model is a system of degenerate reaction-diffusion equations, coupled with a system of ODEs, i.e. with the ionic model of the cell membrane. The derivation of Neu-Krassowska was an important step in understanding the exact influence of the cell scale parameters on the macroscale behaviour, but it was only a formal one. They did not provide a rigorous mathematical proof either for the convergence in their asymptotic analysis, nor for the well-posedness of the bidomain model.

In 2005, Pennacchio, Savare and Colli Franzone [START_REF] Pennacchio | Multiscale modeling for the bioelectric activity of the heart[END_REF], using the framework of the Γconvergence theory and the same assumptions on the microstructure of the cardiac tissue, presented the rigorous mathematical derivation of the bidomain model, as a limit problem of the microscopic (cell scale) model. As for the well-posedness of the bidomain model there are two main challenges: the degeneracy of the reaction-diffusion equations and the non-linearity of the coupled ionic models. The first proof of existence was given in 2002 by Colli Franzone and Savare [START_REF] Colli | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF], for the FitzHugh Nagumo ionic model. The following was the work of Bourgault, Coudière and Pierre in 2009 [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF], where the compactness technique is used and the proof of the existence was extended to the more complex ionic models, namely Aliev-Panfilov and MacCulloch models. The same year Veneroni proved the existence and uniqueness for Luo-Rudy ionic model, [START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF]. In 2008, Boulakia et al. [START_REF] Boulakia | A coupled system of PDEs and ODEs arising in electrocardiograms modeling[END_REF], gave a proof for a coupled heart-torso problem, and extended the proof to include the Mitchell-Schaeffer ionic model.

Bidomain model

Here, we present the intuitive understanding of the bidomain model as it has been described in the 70's. In the microstructure the intracellular and the extracellular spaces exist as separate domains, that interact through the cell membrane. On the macroscale in the averaged model, the tissue is seen as a two-phase medium, as if in every point in space there exist both the intracellular and the extracellular spaces. Hence, for each point in space we can define the intracellular and extracellular potentials, u i and u e , as well as the two currents I i and I e , and finally one can define the transmembrane potential in each point in space, v = u i -u e .

Ohm's law states that the current flux in the conducting body is given with

J = σE, (1.48)
where σ is the conductivity tensor of the medium, and E is the electric field. Under the quasi-static assumption, the electric field E is defined as the gradient of a scalar potential u, i.e. E = -∇u, giving J = -σ∇u.

(1.49)

It holds for both, intracellular and extracellular spaces, respectively,

J i = -σ i ∇u i , ( 1.50
)

J e = -σ e ∇u e , (1.51)
where σ i and σ e are intracellular and extracellular conductivity tensors. The principal axes of the conductivity tensors are the same, but the conductivities in these directions are possibly different. At any point in space the total current flux is J t = J i + J e , and there is no accumulation of charge anywhere, and therefore

∇ • J t = 0, (1.52)
that gives,

∇ • (σ i ∇u i ) + ∇ • (σ e ∇u e ) = 0. (1.53)
The transmembrane current I T , is the current that leaves the intracellular space to enter the extracellular space, so

I T = ∇ • (σ i ∇u i ) = -∇ • (σ e ∇u e ).
(1.54)

For the biological membrane the total transmembrane current, as we mentioned before, is the sum of ionic and capacitive currents,

I T = χ C m ∂v ∂t + I ion , (1.55)
where C m is the specific capacitance per unit area of the cell membrane. Here, χ the surface-to-volume ratio, needed to convert the transmembrane current per unit area into the transmembrane current per unit volume. Boundary conditions for this model usually assume that there is no current across the boundary, i.e. it assumes that the tissue is electrically isolated. The full bidomain model is given as follows: let Ω be a heart tissue (full heart or piece of the tissue) and we want to model the activity in this domain for the time [0, T ], then the bidomain model reads,

∂h ∂t + g(v, h) = 0, in (0, T ) × Ω, χ C m ∂v ∂t + I ion (v, h) = ∇ • (σ i ∇u i ), in (0, T ) × Ω, χ C m ∂v ∂t + I ion (v, h) = -∇ • (σ e ∇u e ), in (0, T ) × Ω, (σ i ∇u i ) • n = (σ e ∇u e ) • n = 0, on [0, T ] × ∂Ω, v(0, x) = v 0 (x), in Ω, h(0, x) = h 0 (x), in Ω, (1.56)
with h the vector of variables for the ionic model. This is the parabolic-parabolic form.

The alternative form of the bidomain model, that is often used in the literature is the parabolic-elliptic form,

∂h ∂t + g(v, h) = 0, in (0, T ) × Ω, χ C m ∂v ∂t + I ion (v, h) = ∇ • (σ i ∇(u e + v)), in (0, T ) × Ω, ∇ • ((σ i + σ e )∇u e ) + ∇ • (σ i ∇v) = 0, in (0, T ) × Ω, (σ i ∇(u e + v)) • n = (σ e ∇u e ) • n = 0, on [0, T ] × ∂Ω, v(0, x) = v 0 (x), in Ω, h(0, x) = h 0 (x), in Ω, (1.57)

Fibers and anisotropy

The conductive properties of the heart muscle are strongly anisotropic. The anisotropy results from the fact that the heart muscle consists of fibers, and the conductivity is higher in the direction of fibers than in the cross-fiber direction. Furthermore, the fibers are organised in sheets, which gives three characteristic directions for the conductive values of the tissue: parallel to the fibers, perpendicular to the fibers and parallel to the sheets, and perpendicular to the sheets. The fiber directions vary through the heart muscle. At a given point we may define a set of unit vectors: a l , a t , and a n , where a l is directed parallel to the fibers, a t is perpendicular to the fibers in the sheet plane, and a n is perpendicular to the sheet plane. The conductivity values in the bidomain model are represented at each point in space by a tensor σ. The local conductivity tensor σ expressed in the eigenbasis of unit vectors a l , a t and a n , is diagonal:

σ =     σ l 0 0 0 σ t 0 0 0 σ n     .
(1.58)

Let A(x) be a local orthogonal matrix having the vectors a l , a t and a n as columns. Then, the conductivity tensor in the global coordinate system is given as

σ(x) = A(x)σA(x) T = σ l a l a T l + σ t a t a T t + σ n a n a T n . (1.59)
Values of the conductivities eigenvalues σ l , σ t and σ n for different media (intracellular, extracellular media, and any other extracardiac media) can be found in the literature (e.g. see Table 3.1). Knowing these values and the proper orientation of the fibers, the conductivity tensor in the global coordinate system is easy to compute. Namely, an entry σ kj in (1.59) is given by:

σ kj (x) = a k l a j l σ l + a k t a j t σ t + a k n a j n σ n , (1.60)
for k, j = 1, 2, 3. Often there is an assumption of the axially symmetric anisotropy, where diffusion in all directions orthogonal to the fibre direction is assumed to be the same, i.e. σ t = σ n . This can be due to the lack of information on the difference in values of σ t and σ n , or because the sheet structure is not visible in the images, or for simplicity of the computations. Using the fact that

a l a T l + a t a T t + a n a T n = 1, (1.61)
and rewriting the equation (1.59) as,

σ = σ n I + (σ l -σ n )a l a T l + (σ t -σ n )a t a T t , (1.62)
the global conductivity tensor is simplified to,

σ = σ n I + (σ l -σ n )a l a T l , (1.63)
or differently,

σ kj = σ n δ kj + (σ l -σ n )a k l a j l , (1.64)
for k, j = 1, 2, 3, where

δ kj = 1 if k = j, otherwise δ kj = 0.
When the conductivity tensors are expressed in their diagonal form, we define the intracellular and the extracellular anisotropy ratios as:

AR int = σ i l σ i n , AR ext = σ e l σ e n .
(1.65)

We will use this definition later in the Chapter 3.

In the ventricles, the fiber orientation has been shown to smoothly rotate from the endocardium to the epicardium, while in other regions the fiber orientation is much more variable. Various techniques can be used to obtain the fiber direction such as histology, optical techniques and diffusion tensor MRI (DT-MRI). An important experimental study on the fiber and sheet structure in ventricular tissue was done by Hooks et al. in 2007, [55], where they demonstrate that the conductivities in the cardiac tissue in fiber direction, parallel, and normal to sheets have the ratio 4:2:1, respectively, which means that the cardiac tissue should be considered as an orthotropic medium. Prior to this study it was widely believed that the the cardiac anisotropy is such that the conductivity is greatest in the fiber direction, and that that it is isotropic in the transverse direction.

Other relevant works in this domain include those of Gilbert et al. [START_REF] Gilbert | Slowed propagation across the compacta-trabeculata interface: a consequence of fiber and sheet anisotropy[END_REF], and Benson et al. [START_REF] Benson | Effects of geometry and architecture on re-entrant scroll wave dynamics in human virtual ventricular tissues[END_REF]. In their studies they use the fiber orientation obtained with DT-MRI of the rat heart, and applied it in the computational model to recover the pattern and velocity of the signal propagation, where they once more confirm the importance of the anisotropic, i.e. orthotropic structure of the tissue.

Asymptotic ionic model

In this subsection we refer to work of M. Rioux and Y. Bourgault, [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF][START_REF] Rioux | Numerical computations of action potentials for the heart-torso coupling problem[END_REF]. The problem at hand is applying the ionic model in the model of the propagation, i.e. in our case into the bidomain model. One needs to match the time and space scales and tune several parameters properly. The idea behind their work is that one can perform the non dimensionalisation and reduce the number of parameters that have to be set to a single nondimensional number. Secondly, using an asymptotic analysis one is able to reproduce the overall AP dynamics using only Mitchell-Schaeffer ionic model with adjusted parameters, based on a knowledge of the resulting solution. Here we give a sketch of their work.

We start from the bidomain model given in (1.57), where we assume it to be coupled with any ionic model of the form,

I ion = h p(v, u e ) τ (1.66)
where h is a gating variable, p(v, u e ) is any type of function of the transmembrane potential and τ is the time scale of this charge flow, in seconds. The first step in the non-dimensionalisation process is to define non-dimensional variables and parameters, denoted by •. The intracellular, extracellular and transmembrane potential are rescaled as:

u i = V m u i + 2V rest , u e = V m u e + V rest and v = V m v + V rest , where V m = V max -V rest ,
and V max and V rest are maximal and resting transmembrane potentials, respectively. V m is called a characteristic transmembrane potential amplitude. The eigenvalues of the conductivity tensors are rescaled as σ i l,t,n = ς σ i l,t,n , and σ e l,t,n = ς σ e l,t,n . The independent variables, time and space, are rescaled as t = T t and x = L x. The gating variables h are assumed to be already non-dimensional as for the most of ionic models in the literature. The non-dimensional version of the bidomain model then reads,

∂h ∂t + g( v, h) = 0, in (0, T ) × Ω, ∂ v ∂ t + I ion ( v, h) = N ∇ • ( σ i ∇( v + u e )), in (0, T ) × Ω, ∇ • ( σ i ∇ v) + ∇ • (( σ i + σ e ) ∇ u e ) = 0, in (0, T ) × Ω, ( σ i ∇( u e + v)) • n = ( σ e ∇ u e ) • n = 0, on [0, T ] × ∂Ω, (1.67)
where,

N = Tς χC m L 2 .
(1.68)

The parameters χ, C m , V m , V max , V rest and ς are fixed. The parameters T and L have to be chosen based on the best suitable time and space scale for the given model and the phase of AP. In the work of Rioux and Bourgault, the suggested parameters are set for the healthy human heart for the MS model, and they are given in the Table 1.1.

Tab. 1.1: Parameters of the bidomain model for a human heart. Three different dimensional scalings are proposed. Adim1: based on the transmembrane potential upstroke. Adim2: based on the action potential duration and width. Adim3: based on convenient units. The table is taken from [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF].

Parameter Adim1 Adim2 Adim3 units

V max 0.04 V V rest -0.085 V V m 0.125 V χ 2 • 10 5 m -1 C m 0.01 F m -2 ς 0.1 S m -1 σ i l 1.741 - σ i t,n 0.1934 - σ e l 3.906 - σ e t,n 1.970 - T 1 • 10 -3 0.2 1 • 10 -3 s L 5 • 10 -5 0.1 1 • 10 -3 m N 0.2 0.001 0.05 -
In the second part of their work they describe the asymptotic analysis performed in order to justify the possibility to use the modified MS model instead of any other ionic model, if we know the time scales of the four phases of the AP signal of interest. More precisely, the parameters that need to be recalculated, for the fixed non-dimensionalisation, i.e. for the fixed N, T and L, are the time constants of the MS model: τ in , τ out , τ open and τ close .

A detailed analysis of the phase plane for the MS model is performed. The phase plane is split into four phases that correspond to depolarisation, AP duration, repolarisation and recovery period. For each phase it is possible to perform the asymptotic analysis and obtain the simplified version of the model.

First, one needs to express the rescaled model in terms of the short and long time scales, i.e. T f and T g , respectively.

∂h ∂t + T T g g( v, h) = 0, in (0, T ) × Ω, ∂ v ∂ t + T T f I ion ( v, h) = N ∇ • ( σ i ∇( v + u e )), in (0, T ) × Ω.
(1.69)

Then, for each phase the rescaling is done again with the time scale and space scale corresponding to the given phase. The time scales are denoted with Tup , TAP , Tdown and Trec . In such settings, some of the terms become negligible and the previous equations are simplified. From these simplified equations and from the properties of the phase plane, one is able to express given time scales with respect to τ 's. The dependencies are the following

Tup = Tup (τ in , τ out ), TAP = TAP (τ close , τ in , τ out ), Tdown = Tdown (τ out ), Trec = Trec (τ open , τ in , τ out ).
(1.70) Each of the above expressions assume fixed v gate , that represents the cell's excitable threshold. Now, since we do have an information on the duration of each phase, we are able to obtain the values for τ 's from the previous dependences in the following order: from Tdown we obtain τ out , then from Tup we obtain τ in . Finally, from TAP and Trec we obtain τ close and τ open , respectively. For all the details on the asymptotic analysis we refer the reader to the papers of Rioux and Bourgault.

With this approach one can replace any ionic model of the type (1.66) with the modified MS model, as far as only the AP overall shape and the time scale are of interest. As said before, the exact ionic models can have dozens of the equations that are hard to solve. This approach in modelling the ionic dynamics can save us a precious computing time. This asymptotic analysis will be of a great use for us in 3D simulations of the rat heart in the Chapter 4.

Limitations of the bidomain model and structural heterogeneities

While the bidomain model is widely accepted as a standard model for the cardiac electrophysiology, it has several limitations, that come from the assumptions in it's derivation. As pointed out in [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF], some of the assumptions can have the important implications, such as: the heart is treated as a static continuum, where parameter values are either uniformly distributed or vary smoothly in space, the cardiac tissue is comprised of myocytes and extracellular space only, while other types of cells and compartments are neglected. Similarly, assumptions on modelling parameter values can have important consequences for the accuracy and validity of models, yet many parameter values are difficult to obtain and verify.

The cardiac tissue is composed of many cell types, supported by the extracellular matrix and permeated by fluids. Myocytes are the most studied cells, they occupy the largest volume of the cardiac tissue, they are electrically active and contractile cells.

On the other hand, the most numerous cells in the heart are fibroblasts. They are much smaller cells than myocytes, and they play an important role in maintenance of the extracellular matrix of the cardiac tissue as they produce interstitial collagen. They are involved in the development of fibrosis in the injured or aged heart. Pathological states are frequently associated with myocardial remodelling involving fibrosis. This is observed in ischemic and rheumatic heart disease, inflammation, hypertrophy, and infarction, [START_REF] Camelliti | Structural and functional characterisation of cardiac fibroblasts[END_REF]. The structural arrangement of the fibroblasts is still not well understood. Some work has been done to understand the role of the fibroblasts in electrophysiology [START_REF] Camelliti | Structural and functional characterisation of cardiac fibroblasts[END_REF][START_REF] Kohl | Mechanosensitive fibroblasts in the sino-atrial node region of rat heart: interaction with cardiomyocytes and possible role[END_REF][START_REF] Kohl | Electrical coupling of fibroblasts and myocytes: relevance for cardiac propagation[END_REF][START_REF] Andrew Maccannell | A mathematical model of electrotonic interactions between ventricular myocytes and fibroblasts[END_REF][START_REF] Frank B Sachse | A model of electrical conduction in cardiac tissue including fibroblasts[END_REF], but it still remains to be studied.

The extracellular matrix is a complex network of fibrous proteins, mainly collagen and elastin. The arrangement of collagen differs throughout the heart. Collagen surrounds each myocyte cell, and it envelops groups of adjacent myocytes and provide the laminar structure of the myocardium. The thickness of collagen fibrils can increase in pathological cases from 40 nm to up to 300 nm, [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF].

There are several open questions in modelling cardiac electrophysiology. They include the choice of the parameters, such as the tissue conductivities, finding a method to represent the detailed and heterogeneous tissue microstructure and the representation of pathological structure and function, [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF].

The gap junctions

The excitation in the working myocardium spreads from the initiating sites to all other cells, providing the synchronisation of all cells. In this sense the heart is said to be a functional syncytium. We know that each cell is surrounded by a plasma membrane, that behaves as a good electrical insulator. One would expect that the electrical current would not be possible to pass from the interior of one cell to the next cell. The solution lies in gap junctions, the structure that allow the diffusion of ions and metabolites by crating a direct pathway between adjacent cells. One of the main reasons for the slow conduction in the AV bundle fibers is caused mainly by diminished numbers of gap junctions between successive cells in the conducting pathways, so that there is great resistance to conduction of excitatory ions from one conducting fiber to the next. Likewise, the rapid transmission of action potentials by Purkinje fibers is believed to be caused by a very high level of permeability of the gap junctions between the successive cells that make up the Purkinje fibers, [START_REF] Hall | Textbook of medical physiology[END_REF].

Gap junctional structure

The cardiac cells are arranged in a brick-like structure. Each cell has a rectangular shape and it is about 100 µm long and 20 µm wide. As we already said, the cells are enclosed in the lipid bilayer -cell membrane, that is about 7 -10 nm wide. The cells are normally at a distance of about 20 nm. Cardiac myocytes are mechanically connected by junctions anchored to the cytoskeleton, known as intercalated disks, and are electrically connected by electrical synapses, known as gap junctions, connecting myocytes mostly end-to-end (longitudinal gap junctions) but also laterally (transverse gap junctions). On the intercalated discs the gap between cells reduces to about 3 nm.

Gap junctions (GJ) are clusters of gap junctional channels (GJCs) that connect the interior of one cell directly to the interior of the neighbouring cells. Gap junction channels are ionic channels with low, mostly size driven selectivity, with a diameter of about 2 nm and length of about 3-15 nm, see Figure 1.10. Each connexon is running transverse to the opposing membrane. Connexons of the adjacent cells align themselves to form the structural continuity. Edited image taken from [START_REF] Makowski | Gap junction structures: Analysis of the x-ray diffraction data[END_REF].

Each GJC is formed by two end-to-end hexameric structures called hemichannels or connexons, each of which is inserted in the membrane of adjacent cells. The hemichannels are composed of connexins, a protein family of 20 isoforms, [START_REF] Harris | Emerging issues of connexin channels: biophysics fills the gap[END_REF]. The cells can express more than one type of connexin, which provides the possibility of forming homomeric and heteromeric connexons, which in turn can give homotypic or heterotypic gap junction channel, shown in the Figure 1.11. Each type of channel exhibit its own set of properties, [START_REF] Desplantez | Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties[END_REF]. Primarily three different connexins: Cx43, Cx40, and Cx45, are expressed in the heart. Cx40 is mainly found in atria and Cx45 in the conduction system. In adult human ventricles, GJCs contain primarily Cx43 [START_REF] Habo | Gap junctions in cardiovascular disease[END_REF], and are located mainly on the longitudinal ends of the cells (Figure 1.12), yielding a macroscopic electrical coupling that is anisotropic [START_REF] Colli Franzone | [END_REF]. Fourteen different connexons can be formed when two connexins are co-expressed in a cell: 2 homomeric connexons, consisting of one type of connexin, and 12 heteromeric connexons consisting of two types of connexins. Gap junction channels formed after docking of two connexons can be homotypic when the connexons are identical (composed of one or two types of connexins) or heterotypic when the connexons are different (composed of one or two types of connexins). Taken from [START_REF] Desplantez | Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties[END_REF]. The anisotropic conductive properties of cardiac tissue are dependent on the geometry of the interconnected cells and the number, size, and location of the gap junction plaques between them, [START_REF] Habo | Gap junctions in cardiovascular disease[END_REF].

Normal adult human ventricular myocardium has a gap junctional surface area of 0.0051 µm 2 /µm 3 myocyte volume, which is comparable to the surfaces in other species rabbit ventricle 0.017 µm 2 /µm 3 , rat ventricle 0.005 -0.015 µm 2 /µm 3 , canine ventricle 0.0085 µm 2 /µm 3 , [START_REF] Peters | Reduced content of connexin43 gap junctions in ventricular myocardium from hypertrophied and ischemic human hearts[END_REF]. There is a correlation between the reduction in this surface area and the conditions of chronic hypertrophy and ischemia, [START_REF] Peters | Reduced content of connexin43 gap junctions in ventricular myocardium from hypertrophied and ischemic human hearts[END_REF].

Ventricular myocytes are connected by gap junctions to about 10 neighbouring cells. GJs are arranged in hexagonal arrays with lattice constants of around 8.0-9.5 nm [START_REF] Wilders | Limitations of the dual voltage clamp method in assaying conductance and kinetics of gap junction channels[END_REF], see Figure 1.12. Conduction velocity is determined by GJ plaque area. Mean gap junctional plaque area ranges from 0.21 µm 2 in the human ventricle, to about 0.45 µm 2 in the rat ventricle, and to 4 µm 2 in the canine ventricle [START_REF] Habo | Gap junctions in cardiovascular disease[END_REF].

Electrical function

In the usual continuous models for the electrical propagation, the GJCs connect cells which interior domains extend the intra-cellular domain [START_REF] Neu | Homogenization of syncytial tissues[END_REF], e.g. simple linear electrical conductors.

However, the dual patch clamp experiments, e.g. [START_REF] Harris | Kinetic properties of a voltage-dependent junctional conductance[END_REF][START_REF] Del Corsso | Transfection of mammalian cells with connexins and measurement of voltage sensitivity of their gap junctions[END_REF][START_REF] Desplantez | Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties[END_REF], have shown that GJCs exhibit a non-linear voltage and time dependence on the difference of membrane potential between adjacent cells. A dual whole cell patch clamp experiment is devised to measure directly the junctional current. Each cell of a coupled pair is impaled with one patch pipette, see Figure 1.13. The voltage steps are delivered to one cell (e.g. cell B), while the other ( e.g. cell A), is held at the constant potential. In such situation the current flowing from the cell B to the cell A, junctional current I j , is equal to the current of the opposite polarity that needs to be injected into cell A to keep it at the constant potential. From here one is able to derive the conductance of the gap junction.

More precisely, let us denote V A and V B the holding potentials for the cells A and B, respectively, I A and I B the currents supplied by voltage clamps, and G A and G B are the conductances of the nonjunctional membrane. Then, the current flowing from the cell A to cell B is

I j = G j V j , where V j = V A -V B .
Then the currents are given as

I A = G A V A + G j (V A -V B ),
and

I B = G B V B + G j (V B -V A ).
When V A is kept at zero and V B is stepped, then I A = -G j V B and V j = -V B , hence the measured current I A = I j , and

G j = I A V j .
Such experiments show that the junctional conductance is time and voltage dependent. For small voltage steps applied, the junctional current and conductance remain almost constant. For the larger steps the junctional current is high in the beginning and decays exponentially over time to the steady state level. It has been observed that the junctional current never reaches zero, meaning there is a residual conductance of gap junctions, G j,min . The initial current and conductance are referred to as instantaneous, and denoted as I j,0 and G j,0 , respectively. The "final" current and conductance are referred to as steady state, and denoted as I j,∞ and G j,∞ . Both, instantaneous and steady state current depend on the voltage step applied.

The dual patch clamp experiment has to be used with care. An interesting study by [START_REF] Wilders | Limitations of the dual voltage clamp method in assaying conductance and kinetics of gap junction channels[END_REF] [START_REF] Wilders | Limitations of the dual voltage clamp method in assaying conductance and kinetics of gap junction channels[END_REF], gives us a deeper insight into this technique and its limitations. They have shown that the actual voltage across each channel of gap junction, in percent of the applied voltage, decreases with the gap junction size. This in return implies that the gap junctional current is not simply proportional to the number of channels in the gap junctions. Hence, the instantaneous gap junctional conductance also, will not be linearly dependent on the number of gap junctional channels, and could be underestimated in such experiments. Additionally, the steady state conductance, G j,∞ , depends on the number of channels in the gap junction. This dependence is much more emphasised if the pipette resistance is high.

Another important observation of this study is that one can not observe the voltage dependence of I j in certain experimental settings. According to this model, there are three parameters that can mask the voltage dependence of the GJCs: the instantaneous conductance might be too large (> 8-10 nS), clamp pulses might be too small (<50 mV) or too short(<1 s), and pipette resistances might be too large (> 10M Ω). This may lead to the wrong conclusion that the gap junction conductance is constant. Finally, this model also gives an explanation for the residual conductance G j,min , as a direct consequence of the non-zero single channel open probability.

Structural heterogeneities

" "If you want to tell someone that you will love them forever you can give them a diamond. But

Introduction

The standard macroscopic model for the electrophysiology of the heart is the bidomain model. It is the anisotropic three-dimensional cable equation, that represents the averaged electric behaviour of the heart tissue. In the bidomain model we have different electrical conductivities for the intracellular and the extracellular spaces, and both of them are anisotropic, meaning that there is a different conductivity in the longitudinal and transversal direction with respect to fiber direction.

Even if it is widely accepted, this model still has several modelling limitations. It assumes that active myocytes are present everywhere, organised into a dense network. It is a reasonable assumption for describing the propagation of the action potential in healthy tissues. As we have discussed in the Section 1.3.4, in pathological cases this assumption does not hold, and there are regions with larger patches of collagen or fibrosis. This is observed in ischemic and rheumatic heart disease, inflammation, hypertrophy, and infarction [START_REF] Camelliti | Structural and functional characterisation of cardiac fibroblasts[END_REF].

In the current modelling of such defects, usually the standard models for the healthy tissues are used, with the model parameters tuned ad hoc. We propose an explanatory model for more rigorous tuning of the parameters in such situations.

In the proposed model of the electrical activity of the heart tissue we will assume the periodic alternations of the healthy tissue, modelled with the standard bidomain model, and the non-active regions, modelled with the simple diffusion equation. Further, we will use the homogenisation technique to derive the macroscale model that can be used in numerical simulations on the whole heart.

Microscopic problem statement 2.2.1 Settings

We consider the bounded open set Ω ∈ R N (N = 2 or 3), with the Lipschitz boundary ∂Ω, such that

Ω = Ω B ε ∪ Ω D ε ∪ Σ ε .
Here Ω B ε represents the healthy heart tissue which can be modelled with standard bidomain equations, Ω D ε represents the collection of periodical diffusive inclusions and Σ ε is the common boundary of these two subdomains. The domain Ω is a periodic medium, i.e. it is divided into the small cells identical to each other. This small cells are identical up to a translation and rescaling by ε to the unit cell Y = [0, 1] N . Furthermore, the unit cell is decomposed in two parts: Y B represents the tissue that can be modelled by the standard bidomain model, Y D is the diffusive inclusion, hence Y = Y B ∪ Y D ∪ Γ where Γ is the interface. Let ε be a sequence of a strictly positive real numbers which tend to zero. We can write

Ω B ε = z∈Z N ε(z + Y B ) ∩ Ω, Ω D ε = z∈Z N ε(z + Y D ) ∩ Ω, Σ ε = z∈Z N ε(z + Γ) ∩ Ω, as on the Figure 2.1. ε Ω Ω B ε Ω D ε Σ ε ∂Ω Y Y D Y B Γ 1

The system of equations for the microscopic problem

The idea is to extend the standard bidomain model on Ω B ε with the periodic diffusive inclusions on Ω D ε and to study their effect on the macroscopic level. Standard bidomain model involves the intracellular electric potential u i ε (t, x) and the extracellular one u e ε (t, x). While the transmembrane potential is denoted as

v ε = u i ε -u e ε .
The bidomain model assumes that the change in the transmembrane voltage results from the ionic activity and the diffusion of the electric potential. In (0, T ) × Ω B ε , holds

∂ t h ε (t, x) + g(v ε (t, x), h ε (t, x)) = 0, (2.1)
∂ t v ε (t, x) + I ion (v ε (t, x), h ε (t, x)) -∇ • σ i ∇u i ε (t, x) = 0, (2.2)
∂ t v ε (t, x) + I ion (v ε (t, x), h ε (t, x)) + ∇ • (σ e ∇u e ε (t, x)) = 0, (2.3)
where σ i and σ e are the time-independent intracellular and extracellular conductivity tensors. They are assumed to be symmetric and positive definite matrices. Their coefficients are assumed to be the periodic functions of the period Y . The ODE in (2.1) with the function I ion represents the ionic model related to the behaviour of the myocardium cells membrane. In our new model we couple this standard bidomain model with the simple diffusion equation in the periodic inclusions:

∇ • σ d ∇u d ε (t, x) = 0, in (0, T ) × Ω D ε . (2.4)
The standard transmission conditions for the bidomain model with the diffusive domain are given on the interface. On [0, T ] × Σ ε , holds

σ i ∇u i ε (t, x) • n Ω B ε = 0, (2.5) (σ e ∇u e ε (t, x)) • n Ω B ε = σ d ∇u d ε (t, x) • n Ω B ε , (2.6) u e ε (t, x) = u d ε (t, x), (2.7)
where n Ω B ε is the unit normal vector from Ω B ε to Ω D ε . On the outer boundary we have homogeneous Neumann conditions on ∂Ω

(σ i ∇u i ε (t, x)) • n = -(σ e ∇u e ε (t, x)) • n = 0, on [0, T ] × ∂Ω. (2.8)
One can notice that the problem is defined up to the same constant for u i ε , u e ε and u d ε . This can be fixed by enforcing the Gauge condition on u e ε , ∀t ∈ [0, T ],

Ω u e ε (t, x)dx = 0.

(2.9)

We assume the initial conditions on v ε and h ε to be

v ε (0, x) = v 0 ε (x) in Ω B ε , h ε (0, x) = h 0 ε (x) in Ω B ε , (2.10)
where

v 0 ε , h 0 ε ∈ H 1 (Ω B ε ).
The similar set of equations have been studied previously in [START_REF] Boulakia | A coupled system of PDEs and ODEs arising in electrocardiograms modeling[END_REF] on the heart -torso problem, where the torso was represented by the diffusive part and the heart, i.e. the bidomain equations were embedded inside of the diffusive domain.

The transmission conditions (2.6) and (2.7) provide the continuity of the potential and the flux of extracellular potential in Ω B ε and the potential in Ω D ε . We can define a new function u ε on the whole domain Ω such that

u ε (t, x) = u e ε (t, x), in Ω B ε , u d ε (t, x), in Ω D ε . (2.11) Note that if u e ε ∈ H 1 (Ω B ε ) and u d ε ∈ H 1 (Ω D ε ) then u ε ∈ H 1 (Ω).
We also define the conductivity tensor σ as

σ(y) = σ e , in Y B , σ d , in Y D .
(2.12)

In this way we can write the equations (2.3) and (2.4) as one, defined on the whole domain Ω and drop the transmission condition (2.7). Redefine

v ε = u i ε -u ε | Ω B ε
and finally, we can summarise the problem to be solved

∂ t h ε + g(v ε , h ε ) = 0, in (0, T ) × Ω B ε , ∂ t v ε + I ion (v ε , h ε ) -∇ • σ i ∇u i ε = 0, in (0, T ) × Ω B ε , 1 Y B (x/ε) (∂ t v ε + I ion (v ε , h ε )) + ∇ • (σ (x/ε) ∇u ε ) = 0, in (0, T ) × Ω, σ i ∇u i ε • n Ω B ε = 0, on [0, T ] × Σ ε , σ e -σ d ∇u ε • n Ω B ε = 0, on [0, T ] × Σ ε , σ i ∇u i ε • n = 0, on [0, T ] × ∂Ω, (σ e ∇u ε ) • n = 0, on [0, T ] × ∂Ω, Ω u ε = 0, in [0, T ], v ε (0, x) = v 0 ε (x) in Ω B ε , h ε (0, x) = h 0 ε (x) in Ω B ε , (2.13)
where the function 1 Y B (y) is the characteristic function of the domain Y B . We have presented several ionic models used to present the electrical activity of one cell in the introduction. For our study we will concentrate on the regularised version of the Mitchell-Schaeffer (MS) model, which is given by the following formulas:

I ion (v, h) = 1 τ in hv 2 (v -1) - 1 τ out v, (2.14) g(v, h) = 1 τ cl + τ cl -τ op τ cl τ op h ∞ (v) (h -h ∞ (v)). (2.15)
The function h ∞ is given by

h ∞ (v) = 1 2 1 -tanh v -v gate η gate , ( 2.16)
and all the parameters are positive constants. MS model has four time constants that correspond to the four phases of the action potential in the cardiac cells.

Formal derivation of the macroscopic problem

As mentioned earlier, the simulation of the microscopic problem (2.13) is very demanding, since we would have to do the space discretisation that is smaller than size of the periodic cell, i.e. ε. For this reason we are rather interested in finding the averaged macroscopic model over the whole domain Ω. This is done by use of the homogenisation technique, [START_REF] Evans | Partial differential equations[END_REF][START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF].

From a mathematical point of view, we have a family of partial differential operators L ε (with coefficients oscillating with period ε), and a family of solutions u ε which, for a given domain Ω and source term f , satisfy,

L ε (u ε ) = f, in Ω (2.17)
complemented by appropriate boundary conditions. Assuming that the sequence u ε converges, in some sense, to a limit u, we look for a so-called homogenized operator L such that u is a solution of

L(u) = f, in Ω (2.18)
Passing from (2.17) to (2.18) is the homogenization process [START_REF] Allaire | Homogenization and two-scale convergence[END_REF]. To find a form for the limit problem, we use ansatz on u ε

u ε (x) = u 0 (x, x/ε) + εu 1 (x, x/ε) + ε 2 u 2 (x, x/ε) + • • • , (2.19)
where each term u k (x, y) is periodic in y. After the substitution, we will be able to identify the cascade problems

L 0 (u 0 ) = f, (2.20) L 1 (u 1 , u 0 ) = 0. (2.21)
In principle, the precise parameters in L 0 can be computed with the use of the second problem, that enables us to express u 1 in the terms of u 0 , and to formulate so called cell problems in the unit cell. In principle, averaging (2.20) with respect to y gives the homogenised problem (2.18). This two-scale asymptotic expansion method is used only to guess the form of the homogenized operator L, as we need a further arguments to prove that the ansatz (2.19) holds. For further details see [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Evans | Partial differential equations[END_REF].

In this section the method of formal two-scale asymptotic expansions is performed in order to find the macroscopic homogenized problem for (2.13). This formal derivation is interesting from the modelling point of view, as it makes it possible to intuitively obtain the homogenized problem. Assume that the solution (u i ε , u ε , h ε ) can be written as a series

u i ε (t, x) = ∞ k=0 ε k u i k (t, x, x/ε), u ε (t, x) = ∞ k=0 ε k u k (t, x, x/ε), h ε (t, x) = ∞ k=0 ε k h k (t, x, x/ε), (2.22)
with u i k (t, x, y), u k (t, x, y), and h k (t, x, y), Y -periodic functions with respect to the variable y ∈ Y . Using this ansatz, and plugging it into the system of equations (2.13) we obtain the cascade systems of equations with respect to the exponent of ε. We are interested in the first three terms, for k = 0, 1 and 2, that we will use to obtain the macroscopic problem. This is summed up in the following lemma.

Formal lemma 1. Under assumption (2.22), the zero order term in the expansion is the solution of the homogenised problem

∇ x • σ i * ∇ x u i 0 = (∂ t v 0 + I ion (v 0 , h 0 ))|Y B |, in [0, T ] × Ω, ∇ x • (σ * ∇ x u 0 ) = -(∂ t v 0 + I ion (v 0 , h 0 ))|Y B |, in [0, T ] × Ω, ∂ t h 0 + g(v 0 , h 0 ) = 0, in [0, T ] × Ω, ∇u i 0 • n = ∇u 0 • n = 0, on [0, T ] × ∂Ω, v 0 (0, x) = v 0 0 (x), in Ω, h 0 (0, x) = h 0 0 (x), in Ω, (2.23)
and the homogenised (the effective) conductivity tensors are given by its entries, for k, j = 1..N,

σ i * kj = σ i kj |Y B | + σ i k1 Y B ∂ y1 w i j dy + σ i k2 Y B ∂ y2 w i j dy + σ i k3 Y B ∂ y3 w i j dy , σ * kj = Y σ kj dy + Y σ k1 ∂ y1 w j dy + Y σ k2 ∂ y2 w j dy + Y σ k3 ∂ y3 w j dy , (2.24)
and for j = 1..N, w i j , w j are the solutions of the two cell problems, for the intracellular space

∇ • (σ i ∇w i j ) = 0, in Y B , σ i (∇w i j + e j ) • n = 0, on Γ, w i j is Y periodic, (2.25)
and for the extracellular space

∇ • (σ∇w j ) = 0, in Y, (σ e -σ d )(∇w j + e j ) • n = 0, on Γ, w j is Y periodic. (2.26)
Furthermore, u i 1 and u 1 are given as

u i 1 (t, x, y) = N j=1 w i j (y) ∂u i 0 ∂x j (t, x), u 1 (t, x, y) = N j=1 w j (y) ∂u 0 ∂x j (t, x).
(2.27)

Proof. The series (2.22) are plugged into the problem (2.13), and the following derivation rule is used:

∇u k x, x ε = 1 ε ∇ y u k + ∇ x u k x, x ε , ( 2.28)
where ∇ x and ∇ y denote the partial derivative with respect to the first and the second variable of u k (x, y). The same rule holds for u i k and h k . Then the system of equations become a series in ε. On Ω × Y B we have

∂ t h 0 + g(v 0 , h 0 ) + O(ε) = 0,
and

∂ t v 0 + I ion (v 0 , h 0 ) -ε -2 ∇ y • (σ i ∇ y u i 0 ) + ε -1 ∇ x • (σ i ∇ y u i 0 ) + ∇ y • (σ i ∇ x u i 0 ) + ∇ y • (σ i ∇ y u i 1 ) +∇ x • (σ i ∇ x u i 0 ) + ∇ x • (σ i ∇ y u i 1 ) + ∇ y • (σ i ∇ x u i 1 ) + ∇ y • (σ i ∇ y u i 2 ) + O(ε) = 0,
and on Ω × Y 1 Y B (y)(∂ t v 0 + I ion (v 0 , h 0 )) -ε -2 ∇ y • (σ∇ y u 0 ) + ε -1 (∇ x • (σ∇ y u 0 ) + ∇ y • (σ∇ x u 0 ) + ∇ y • (σ∇ y u 1 )) +∇ x • (σ∇ x u 0 ) + ∇ x • (σ∇ y u 1 ) + ∇ y • (σ∇ x u 1 ) + ∇ y • (σ∇ y u 2 )] + O(ε) = 0,
The boundary conditions similarly become series in ε. On Σ ε × Γ, we have

σ i [ε -1 ∇ y u i 0 + ∇ x u i 0 + ∇ y u i 1 + ε(∇ x u i 1 + ∇ y u i 2 )] • n Γ + O(ε 2 ) = 0, (σ e -σ d )[ε -1 ∇ y u 0 + ∇ x u 0 + ∇ y u 1 + ε(∇ x u 1 + ∇ y u 2 )] • n Γ + O(ε 2 ) = 0,
with the periodic boundary conditions on ∂Y . Identifying each coefficient of ε as an individual equation yields a cascade of systems of equations. The ε -2 system of equations gives, for u i 0 ,

∇ y • σ i ∇ y u i 0 = 0, in Ω × Y B , σ i ∇ y u i 0 • n Γ = 0, on Σ ε × Γ,
+ periodic boundary conditions on ∂Y.

(2.29)

If we multiply (2.29) with u i 0 and integrate over y, we obtain

0 = Y B u i 0 ∇ y • σ i ∇ y u i 0 dy = - Y B σ i ∇ y u i 0 • ∇ y u i 0 dy + Γ u i 0 σ i ∇ y u i 0 • n Γ dy. =0 (2.30)
Hence,

Y B |∇ y u i 0 | 2 dy = 0 =⇒ ∇ y u i 0 = 0 =⇒ u i 0 = u i 0 (t, x). (2.31)
And for u 0 we get,

∇ y • (σ∇ y u 0 ) = 0, in Ω × Y, (σ e -σ d )∇ y u 0 • n Γ = 0, on Σ ε × Γ,
+ periodic boundary condition on ∂Y.

(2.32)

Similarly, if we multiply (2.32) with u 0 and integrate over y, we have

0 = Y u 0 ∇ y • (σ∇ y u 0 )dy (2.33) = - Y (σ∇ y u 0 ) • ∇ y u 0 dy + Γ u 0 ((σ e -σ d )∇ y u 0 ) • n Γ dy. =0 (2.34)
Hence,

Y |∇ y u 0 | 2 dy = 0 =⇒ ∇ y u 0 = 0 =⇒ u 0 = u 0 (t, x). (2.35)
Meaning, the terms u i 0 and u 0 do not depend on a small scale y. These terms represent the limit values of u i ε and u ε , i.e. the averaged intra and extracellular potentials. The ε -1 system of equation then gives,

∇ y • (σ i ∇ y u i 1 ) = 0, in Ω × Y B , σ i ∇ y u i 1 • n Γ = -σ i ∇ x u i 0 • n Γ , on Σ ε × Γ, + periodic boundary condition on ∂Y, (2.36)
and,

∇ y • (σ∇ y u 1 ) = 0, in Ω × Y, (σ e -σ d )∇ y u 1 • n Γ = -(σ e -σ d )∇ x u 0 • n Γ , on Σ ε × Γ,
+ periodic boundary condition on ∂Y.

(2.37)

From (2.36) and (2.37) we see that the terms u i 1 and u 1 can be expressed as functions of u i 0 and u 0 . Finally, the ε 0 system of equation gives the equation for the term

u i 2 in Ω × Y B , -∇ y • (σ i ∇ y u i 2 ) =∇ y • (σ i ∇ x u i 1 ) + ∇ x • (σ i ∇ x u i 0 ) + ∇ x • (σ i ∇ y u i 1 ) -∂ t v 0 + I ion (v 0 , h 0 ), (2.38)
with the boundary condition,

σ i ∇ y u i 2 • n Γ = σ i ∇ x u i 1 • n Γ , on Σ ε × Γ (2.39)
+ periodic boundary condition on ∂Y.

(2.40)

Integrating the left hand side of (2.38) over Y B , and using the periodic boundary condition on u i 2 , we obtain

- Y B ∇ y • (σ i ∇ y u i 2 )dy = Γ (σ i ∇ y u i 2 ) • n Γ ds (2.41) = Γ (σ i ∇ x u i 1 ) • n Γ ds (2.42) = - Y B ∇ y • (σ i ∇ x u i 1 )dy, (2.43)
which implies that the right hand side of (2.38), when integrated over Y B , gives

Y B ∇ x • (σ i ∇ x u i 0 ) + ∇ x • (σ i ∇ y u i 1 ) dy - Y B (∂ t v 0 + I ion (v 0 , h 0 )) dy = 0, (2.44)
and thus

|Y B | (∂ t v 0 + I ion (v 0 , h 0 )) = |Y B |∇ x • (σ i ∇ x u i 0 ) + Y B ∇ x • (σ i ∇ y u i 1 )dy. (2.45)
And for the term u 2 in Ω × Y , the system of equation gives

-∇ y • (σ∇ y u 2 ) =∇ y • (σ∇ x u 1 ) + ∇ x • (σ∇ x u 0 ) + ∇ x • (σ∇ y u 1 ) -1 Y B (y) (∂ t v 0 + I ion (v 0 , h 0 )) , (2.46)
with the boundary conditions

(σ e -σ d )∇ y u 2 • n Γ = (σ e -σ d )∇ x u 1 • n Γ , on Σ ε × Γ (2.47)
+ periodic boundary conditions on ∂Y.

(2.48)

Integrating the left hand side of (2.46) over Y , and using the periodic boundary condition on u 2 , we obtain

- Y ∇ y • (σ∇ y u 2 )dy = - Y B ∇ y • (σ e ∇ y u 2 )dy - Y D ∇ y • (σ d ∇ y u 2 )dy (2.49) = Γ ((σ e -σ d )∇ y u 2 ) • n Γ ds (2.50) = Γ ((σ e -σ d )∇ x u 1 ) • n Γ ds (2.51) = - Y ∇ y • (σ∇ x u 1 )dy. (2.52)
Therefore the right hand side of (2.46), when integrated over Y B , gives

Y (∇ x • (σ∇ x u 0 ) + ∇ x • (σ∇ y u 1 )) dy - Y 1 Y B (y) (∂ t v 0 + I ion (v 0 , h 0 )) dy = 0, (2.53)
which simplifies into

|Y B | (∂ t v 0 + I ion (v 0 , h 0 )) = Y ∇ x • (σ∇ x u 0 ) + Y ∇ x • (σ∇ y u 1 )
∂ t h 0 + g(v 0 , h 0 ) = 0,
represent the homogenised macroscopic problem, whose solution (h 0 , u i 0 , u 0 ) is the limit of (h ε , u i ε , u ε ). The small scale effects are accounted for through the terms u i 1 and u 1 , that can be expressed in terms of u i 0 and u 0 , from (2.36) and (2.37). In fact, these terms are defined up to addition of functions that depend only on x, but this does not matter since we only use their gradients ∇ y u i 1 and ∇ y u 1 to obtain the homogenised equations for u i 0 and u 0 . In order to compute u i 1 and u 1 and to simplify the homogenised problem, we introduce the so-called cell problems. We express u i 1 (t, x, y) and u 1 (t, x, y) as in (2.27), and use the equations (2.36) and (2.37) to define the cell problems for functions w i j (y) and w j (y). The simple calculation leads us to the problems in the periodic unit cell given in (2.25) and (2.26). The functions w i j (y) and w j (y) are the local variations of the intacellular an extracellular potentials, created by an averaged gradient e j .

Finally to simplify the homogenised problem, we substitute the expressions (2.27), into (2.45) and (2.54). The resulting homogenised problem is given in (2.23), with the homogenised conductivities given in (2.24). The constant tensors σ i * and σ * describe the effective or homogenised properties of the heterogeneous tissue. They do depend neither on time, nor on the choice of Ω. They depend on the initial conductivity tensors and on the shape and size of the diffusive inclusions.

Remark 2. The method of the two scale convergence, and hence our derivation of the homogenised problem is not a rigorous from a mathematical point of view. In other words, it leads heuristically the homogenised system of equations, but it is not a correct and sufficient proof of the homogenised process. First issue is that we still did not prove the existence of the solution to the original ε-dependent microscopic problem. Secondly, we do not know if such microscopic problem has in fact a limit. And even if we assume it has, the ansatz (2.27) is not yet justified. Hence, we can not claim that (h 0 , u i 0 , u 0 ) is actually a limit to (h ε , u i ε , u ε ) when ε → 0. And the last, but not least issue is that we have to deal with the non-linearity and actually justify the assumption that g(v ε , h ε ) → g(v 0 , h 0 ), and

I ion (v ε , h ε ) → I ion (v 0 , h 0 ).
The details on the rigorous proof are laid out in the Section 2.4.

Numerical convergence study

We want to numerically confirm the convergence of the microscopic problem to the derived homogenised equations, and to observe the rate of the convergence. In order to do this we run simulations of the microscopic problem (2.13), with reducing size of ε, and for the homogenised problem (2.23). To be sure that we observe the change with respect to ε, we run all simulations with the same mesh step, time step, initial conductivities and the Mitchell-Schaeffer mode parameters. To compute the rate of the convergence, we obtain the numerical solutions of the microscopic problems for the time T = 3.0, and compute the L 2 -errors of the microscopic solutions with respect to the solution of the homogenised problem at the same time.

Microscale problem

To run the simulation for the microscopic problem, we create a mesh of the domain with diffusive inclusions, as in 2.1. The u i ε , v ε and h ε are defined on Ω B ε , while u ε is defined everywhere on Ω. For the convergence study we will use only the values of v ε and h ε . ) and u e εn , φ e ∈ H 1 (Ω), for n = 0, .., N. 5: Define v ε0 := v 0 , h ε0 := h 0 . 6: Solve the coupled system with EE to obtain u i ε1 , u e ε1 :

Algorithm 1 Microscale problem

Ω B ε (u i ε1 -u e ε1 )(φ i -φ e ) + ∆t Ω B ε σ i ∇u i ε1 ∇φ i + Ω (σ∇u e ε1 ) ∇φ e + Ω ηu e ε1 φ e = Ω B ε (v ε0 -∆tI ion (v ε0 , h ε0 )) (φ i -φ e ).
7: Solve the ODE to obtain h ε1 :

Ω h ε1 φ = Ω (h ε0 + ∆tg(v ε0 , h ε0 )) φ. 8: Define v ε1 := u i ε1 -u e ε1 | Ω B ε . 9: for n = 1 to N do 10:
Solve the coupled system with SBDF2 to obtain u i εn+1 , u e εn+1 :

Ω B ε 3 2 u i εn+1 -u e εn+1 (φ i -φ e )+∆t Ω B ε σ i ∇u i εn+1 ∇φ i + Ω σ∇u e εn+1 ∇φ e + Ω ηu e εn+1 φ e = Ω B ε (2v εn - 1 2 v εn-1 ) + ∆t (2I ion (v εn , h εn ) -I ion (v εn-1 , h εn-1 )) (φ i -φ e ) .
11:

Solve the ODE to obtain h εn+1 :

Ω B ε 3 2 h εn+1 φ = Ω 2h εn - 1 2 h εn-1 + ∆t (2g(v εn , h εn ) -g(v εn-1 , h εn-1 )) φ.
12:

Define

v εn+1 := u i εn+1 -u e εn+1 | Ω B ε .
Note that here σ is defined as in (2.12), i.e. it takes value σ e on Ω B ε , and σ d on Ω D ε . The term 0 < η 1 is used to impose the Gauge condition on u e ε , given in (2.9).

Homogenised problem

To run the simulation for the macroscopic problem, we need to solve the cells problems first, on the unit cell Y = Y B ∪ Y D , and to compute the modified conductivities. Then we use these to solve the modified bidomain model on the whole domain Ω, without distinguishing Ω B ε and Ω D ε , i.e. without including the boundary Γ ε . The cell problems are static problems, that we solve only once on the beginning of the simulation to obtain the modified conductivities σ i * and σ * .

Algorithm 2 Static cell problems

Y B (σ i ∇w i k ) • ∇φ i + Γ (σ i e k ) • n Γ φ i = 0 5: Solve for w k : Y B (σ∇w k ) • ∇φ + Γ ((σ e -σ d )e k ) • n Γ φ = 0
6: Compute derivatives: ∂w i k and ∂w k . 7: Compute the modfied conductivities σ i * and σ * using (2.24). Now when we have computed the modified conductivities, we can run the simulation for the homogenised problem, (2.23). For simplicity, we omit the subscript 0.

Algorithm 3 Homogenised problem 1: Define mesh on Ω. 2: Define T and N = T ∆t . 3: Define var. form. space H 1 (Ω) with P1 elements. 4: Declare u i n , u n , v n , h n , φ i , φ e , φ ∈ H 1 (Ω), for n = 0, .., N. 5: Define v 0 := v 0 , h 0 := h 0 . 6: Solve the coupled system with EE to obtain u i ε1 , u e ε1 :

|Y B | Ω (u i 1 -u 1 )(φ i -φ e ) + ∆t Ω σ i * ∇u i 1 ∇φ i + Ω (σ * ∇u 1 ) ∇φ e + Ω ηu 1 φ e = |Y B | Ω (v 0 -∆tI ion (v 0 , h 0 )) (φ i -φ e ).
7: Solve the ODE to obtain h 1 :

Ω h 1 φ = Ω (h 0 + ∆tg(v 0 , v 0 )) φ. 8: Define v 1 := u i 1 -u 1 . 9: for n = 1 to N do 10:
Solve the coupled system with SBDF2 to obtain u i n+1 , u n+1 :

|Y B | Ω 3 2 u i n+1 -u n+1 (φ i -φ e )+∆t Ω σ i * ∇u i n+1 ∇φ i + Ω (σ * ∇u n+1 ) ∇φ e + Ω ηu n+1 φ e = |Y B | Ω (2v n - 1 2 v n-1 ) + ∆t (2I ion (v n , h n ) -I ion (v n-1 , h n-1 )) (φ i -φ e ) .
11:

Solve the ODE to obtain h εn+1 :

Ω 3 2 h n+1 φ = Ω 2h n - 1 2 h n-1 + ∆t (2g(v n , h n ) -g(v n-1 , h n-1 )) φ.
12:

Define v n+1 := u i n+1 -u n+1 .

Convergence

To test the convergence numerically, we run both, microscale and macroscale simulation in 2D. We use the Michel Schaeffer ionic model, with the parameters given in the Table 1.1. We choose the domain Ω in the shape of square with the non-dimensional size 60. All simulations are ran on the same mesh steps, dx = 0.3. The final time is T = 3 and the time step dt = 0.5. Further, we assume that the diffusive inclusions are of the circular shape and that they occupy 20% of the total volume, i.e. |Y B | = 0.8. Microscopic problem we run for several sizes of periodic cell, ε ∈ {1/10, 1/15, 1/20, 1/25, 1/30, 1/35, 1/40} . In order to observe numerically the theoretical convergence, with respect to ε, we run each simulation with the same mesh step. We save v ε and h ε , defined on Ω B ε , at the final time T = 10.0. Macroscopic problem is computed using the algorithms 2 and 3. We use the same mesh step as for the microscopic problems, and save the v and h, defined on Ω at the final time T = 10.0. To find the rate of convergence, for each ε we compute the L 2 norm of the differences v ε -v L 2 and h ε -h L 2 , then we compute the norm of v L 2 and h L 2 , and compute so called

L 2 -errors L2err(v ε ) = v ε -v L 2 v L 2 ,
and

L2err(h ε ) = h ε -h L 2 h L 2 .
We use the log-log scale to fit the results with the linear functions. We have obtained approximately linear rates of the convergence. More precisely, we have the rate of convergence 1.39 for v ε , and 0.63 for h ε , as shown in the Figure 2.3.

Fig. 2.3:

Convergence study for v ε and h ε . Observed convergence rate for the given time: 1.39 for V , and 0.63 for h.

The rigorous mathematical proofs: existence and convergence

As we have mentioned earlier, we did the model setting and the convergence study without the rigorous mathematical justification. Here will give all the necessary details on this part.

Existence of the weak solution for the microscopic problem

We go back to our microscopic problem, and we want to prove the existence of the solution for (2.13), by use of the semigroup theory approach [START_REF] Cazenave | An introduction to semilinear evolution equations[END_REF]. In this part we will use the parabolic-elliptic form of the problem.

∂ t h ε + g(v ε , h ε ) = 0, in Ω B ε , (2.55a) ∂ t v ε + I ion (v ε , h ε (t, v ε )) -∇ • (σ i ∇v ε ) = ∇ • (σ i ∇u ε ), in Ω B ε , (2.55b)
-∇ • ((σ i + σ e )∇u ε ) = ∇ • (σ i ∇v ε ), in Ω B ε , (2.55c)
-∇ • (σ d ∇u ε ) = 0, in Ω D ε , (2.55d)
where u ε is continuous. The boundary and transmission conditions are given as

(σ i ∇v ε ) • n Ω B ε = -(σ i ∇u ε ) • n Ω B ε , on Σ ε , (2.55e) (σ e ∇u ε ) • n Ω B ε = -(σ d ∇u ε ) • n Ω D ε , on Σ ε , (2.55f) (σ i ∇v ε ) • n = 0, on ∂Ω, ( 2.55g)
(σ i + σ e )∇u ε • n = 0, on ∂Ω, (2.55h) (σ d ∇u ε ) • n = 0, on ∂Ω, ( 2.55i)
the initial conditions are

h ε (0, x) = h 0 ε (x), in Ω B ε , (2.55j) v ε (0, x) = v 0 ε (x), in Ω B ε , (2.55k)
and the Gauge condition is

Ω u ε dx = 0, (2.55l)
where n Ω B ε is the normal on the Σ ε outwards from Ω B ε , n Ω D ε is the normal on the Σ ε outwards from Ω D ε , and n is the outwards normal on the ∂Ω.

Remark 3. Given the initial conditions for h ε and v ε , and assuming that there exists a unique solution to the above problem (2.55), continuous in time, then u 0 ε (x) := u ε (0, x) is the solution to

Ω B ε (σ i + σ e )∇u 0 ε ∇φ dx + Ω D ε σ d ∇u 0 ε ∇φ dx = - Ω B ε σ i ∇v 0 ε ∇φ dx, ( 2.56)
for all φ ∈ H 1 (Ω), with the Gauge condition Ω u 0 ε dx = 0.

Normally we work with the constant conductivity tensors σ i , σ e and σ d . But our proof will work in the more general case, under the following assumptions on the conductivities.

Assumption 4. We suppose that the conductivity tensors σ e (x), σ i (x) and σ d (x) are symmetric definite positive matrix-functions, such that

σ|ξ| 2 ≤ ξ T σ e,i (x)ξ ≤ σ|ξ| 2 , ∀x ∈ Ω B ε ,
and

σ|ξ| 2 ≤ ξ T σ d (x)ξ ≤ σ|ξ| 2 , ∀x ∈ Ω D ε .
We denote by Assumption 5. We assume that the ionic function I ion (v, h) is globally Lipschitz with respect to v, in a sense that there exist a K > 0, such that,

√ σ i , √ σ e ,
|I ion (v 1 , h(t, v 1 )) -I ion (v 2 , h(t, v 2 ))| ≤ K|v 1 -v 2 |, ∀t ∈ [0, ∞), ∀v 1 , v 2 ∈ R 2 . Lemma 6 (A priori estimates on h ε .). Let h 0 ε ∈ L ∞ (Ω B ε ), with 0 < h 0 ε ≤ 1, then we have h ε (t, •) ∈ L ∞ (Ω B ε ) and ∂ t h ε (t, •) ∈ L ∞ (Ω B ε ). Moreover 0 < h ε (t, x) ≤ 1.
Proof. From the model (2.55) we have

∂ t h ε = -g(v ε , h ε ),
where the function g(v, h) is given in (2.15), so

-g(v ε , h ε ) = (h ∞ (v ε ) -h ε ) 1 τ cl + τ cl -τ op τ cl τ op h ∞ (v ε ) ( * )
.

(2.57)

The term ( * ) is positive, because 0 < τ op < τ cl . From (2.16) we have that

0 ≤ h ∞ (v) ≤ 1, ∀v ∈ R. (2.58)
Hence,

∂ t h ε (t, x) ≥ -h ε (t, x) 1 τ cl + τ cl -τ op τ cl τ op h ∞ (v ε ) , (2.59) ∂ t h ε (t, x) ≤ (1 -h ε (t, x)) 1 τ cl + τ cl -τ op τ cl τ op h ∞ (v ε ) . (2.60)
Then, applying the Gromwall lemma, we have the upper and lower bound on

h ε h ε (t, x) ≥ h 0 ε (x)e - T 0 1 τ cl + τ cl -τop τ cl τop h∞(vε) ds , ( 2.61
)

h ε (t, x) ≤ 1 -(1 -h 0 ε (x))e - T 0 1 τ cl + τ cl -τop τ cl τop h∞(vε) ds . ( 2
Finally, combining (2.61),(2.62) and (2.64), and using the assumption that the initial condition h 0 ε (x) is bounded, we obtain

h min ≤ h ε (t, x) ≤ 1, in Ω B ε , ( 2.65)
where h min := min{h 0 ε (x)}e -T τop . Finally, from (2.59) -(2.60), we have

|∂ t h ε (t, x)| ≤ 1 τ op , in Ω B ε .
(2.66)

In this section, the parameter ε is fixed, so we will omit it as a subscript in the equations in the rest of the section. To build the proof we will use the semigroup approach. Definition 7. We define the operator B is the operator from

H 1 (Ω B ε ) to H 1 (Ω), such that B(v) := u,
and u is the solution of the equation

Ω B ε (σ i + σ e )∇u∇φ dx + Ω D ε σ d ∇u∇φ dx = - Ω B ε σ i ∇v∇φ dx, ∀φ ∈ H 1 (Ω), Ω φ = 0 .
(2.67)

Additionally, we denote as B Ω B ε the restriction of the operator B to Ω B ε :

B Ω B ε (v) := u| Ω B ε .
Lemma 8. For a constant β > 0, the following inequality holds

∇B(v) L 2 (Ω B ε ) ≤ β ∇v L 2 (Ω B ε ) .
(2.68)

Proof. Taking φ = u in (2.67), we have

Ω B ε | √ σ i ∇u| 2 dx + Ω B ε | √ σ e ∇u| 2 dx + Ω D ε | √ σ d ∇u| 2 dx = - Ω B ε σ i ∇v∇udx.
From the Assumption 4, we have

σ ∇u 2 L 2 (Ω B ε ) ≤ σ ∇v L 2 (Ω B ε ) ∇u L 2 (Ω B ε ) .
Taking β = σ/σ, we obtain the inequality (2.68).

Definition 9.

Let A be the operator from

H 1 (Ω B ε ) to H 1 (Ω B ε ) defined as A(v) := -∇ • σ i ∇v -∇ • σ i ∇B Ω B ε (v) , (2.69)
with the domain D(A), Proof. See Theorem 3.4.4. in [START_REF] Cazenave | An introduction to semilinear evolution equations[END_REF].

D(A) := v ∈ H 1 (Ω B ε ) : B Ω B ε (v) ∈ L 2 (Ω B ε ), σ i ∇v + σ i ∇B Ω B ε (v) • n Ω B ε = 0 . ( 2

Lemma 12. The operator (A, D(A))

, is m-dissipative with dense domain in L 2 (Ω B ε ). Therefore the operator A generates a contraction semi-group, whose generator is denoted by e -tA . Proof. To prove that the operator A is m-dissipative, we need to prove that for any λ > 0, and any f ∈ L 2 (Ω B ε ), there exists a unique solution v ∈ D(A) to the equation

v -λA(v) = f, or expanded, v + λ ∇ • σ i ∇v + ∇ • σ i ∇B Ω B ε (v) = f.
The variational formulation for this problem reads, for any

φ B ∈ H 1 (Ω B ε ), Ω B ε vφ B dx + λ Ω B ε σ i ∇v∇φ B dx + Ω B ε σ i ∇B Ω B ε (v)∇φ B dx =a(v,φ B ) = Ω B ε f φ B dx, (2.71)
The boundedness of the bilinear form a(v, φ B ) is shown by using (2.68),

|a(v, φ B )| ≤ v L 2 (Ω B ε ) φ B L 2 (Ω B ε ) + λσ(1 + β) ∇v L 2 (Ω B ε ) ∇φ B L 2 (Ω B ε ) ≤ v H 1 (Ω B ε ) φ B H 1 (Ω B ε ) .
(2.72)

To prove coerciveness we use the definition of the operator B(v) (2.67), and express the bilinear form as follows,

a(v, φ B ) = Ω B ε vφ B dx + λ Ω B ε σ i ∇v∇φ B dx + Ω B ε σ i ∇u∇φ B dx + Ω B ε (σ i + σ e )∇u∇φ dx + Ω B ε σ i ∇v∇φ dx + Ω D ε σ d ∇u∇φ dx .
Then, for φ B = v, and φ = u,

a(v, v) = Ω B ε v 2 dx + λ Ω B ε σ i ∇v∇v dx + Ω B ε σ i ∇u∇v dx + Ω B ε (σ i + σ e )∇u∇u dx + Ω B ε σ i ∇v∇u dx + Ω D ε σ d ∇u∇u dx
then using the inequality: 0 ≤ |a + δb| 2 , ∀a, b, we have

≥ Ω B ε v 2 dx + λ Ω B ε √ σ i ∇v 2 dx -δ Ω B ε √ σ i ∇v 2 dx - 1 δ Ω B ε √ σ i ∇u 2 dx + Ω B ε √ σ i ∇u 2 dx + Ω B ε √ σ e ∇u 2 dx + Ω D ε √ σ d ∇u 2 dx
then using the Assumption 4, and choosing 1 1+σ/σ < δ < 1 we have

≥ Ω B ε v 2 dx + λ Ω B ε (1 -δ)σ |∇v| 2 dx + Ω B ε 1 - 1 δ σ |∇u| 2 + Ω B ε σ |∇u| 2 dx + Ω D ε σ |∇u| 2 dx ≥ Ω B ε v 2 dx + λσ Ω B ε (1 -δ) |∇v| 2 dx + Ω B ε 1 + 1 - 1 δ σ σ |∇u| 2 + Ω D ε |∇u| 2 dx
Now, thanks to the choice of δ we have that every term above is positive. Hence,

a(v, v) ≥ c v H 1 (Ω B ε ) . (2.74)
From Lax-Milgram theorem, the problem (2.71) has a unique solution and both conditions in the Definition 10 are satisfied. Hence, the operator A(v) is m-dissipative, and generates the semigroup denoted by e -tA .

Lemma 13. The evolution problem (2.55) is equivalent to the problem

∂ t h ε + g(v ε , h ε ) = 0, in Ω B ε , ∀t ∈ [0, T ], (2.75a
)

∂ t v ε + A(v ε ) + I ion (v ε , h ε (t, v ε )) = 0, in Ω B ε , ∀t ∈ [0, T ], (2.75b)
with the initial conditions

h ε (0, x) = h 0 ε (x), v ε (0, x) = v 0 ε (x). (2.75c)
Moreover, v ε is also the solution of the following problem:

v ε (t, x) = e -tA v 0 ε - t 0 e -(t-s)A I ion (s, h ε (s, v ε (s, x))) ds, ∀t ∈ [0, T ]. (2.76)
Proof. The equivalence of the two problems comes from the definition of the operators A and B. The formula (2.76) is a straightforward use of the standard Duhamel's formula, see Lemma 4.1.1 in [START_REF] Cazenave | An introduction to semilinear evolution equations[END_REF].

Theorem 14 (Well-posedness of the model). Let ε > 0 be given. Let

h 0 ε ∈ L ∞ (Ω B ε ), with 0 < h 0 ε (x) ≤ 1 and let v 0 ε ∈ L 2 (Ω B ε ). Then the problem 2.55 has a unique solution (v ε , u ε , h ε ), such that v ε ∈ C([0, +∞); H 1 (Ω B ε )) ∩ C 1 ([0, +∞); L 2 (Ω B ε )), u ε ∈ C([0, +∞); H 1 (Ω)), h ε ∈ C([0, +∞); L ∞ (Ω B ε )) ∩ C 1 ([0, +∞); L ∞ (Ω B ε )).
Proof. The proof that v ε exists and is unique in

C([0, +∞); H 1 (Ω B ε ))∩C 1 ([0, +∞); L 2 (Ω B ε )
) is a straightforward consequence of the Hille-Yosida theorem, since (A, D(A)) is a mdissipative operator with dense domain, and the Lipschitz assumption 5 on the ionic function. See the Section 4.3. in [START_REF] Cazenave | An introduction to semilinear evolution equations[END_REF]. The regularity on h ε comes directly from the assumption on the initial conditions and Lemma 6, and the regularity on u ε from the Definition 7 and the Poincare's inequality.

Remark 15. For the function u

i ε := v ε + u ε 1 Ω B ε , we have u i ε ∈ C([0, +∞); H 1 (Ω B ε )).
Remark 16 (On the Lipschitz property of the ionic function). The standard Mitchell-Schaeffer model [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF], which is commonly used for the applications, reads as follows:

∂h ∂t = g(v, h) :=    1-h τopen if v ≤ v gate -h τ close if v > v gate , I ion (v, h(t, v)) = 1 τ in h(t, v)v 2 (v -1) + v τ out .
Thus model has been derived to account for the 4 dynamics of the transmembrane voltage evolution, thus for physiological values of v. Similarly to Boulakia et al. in [START_REF] Boulakia | A coupled system of PDEs and ODEs arising in electrocardiograms modeling[END_REF] we propose here a regularised version of such current, by regularizing the function h. The regularised ionic model Boulakia et al. is given in (2.14) -(2.15), with the regularisation on the step function introduced through h ∞ (v). However, the current remains unbounded, since it grows cubically with respect to v.

One can notice that in the MS model, the quantity

S(v, h(t, v)) = 1 τ in hv(v -1) + 1 τ out ,
is homogeneous to the ratio of the surface conductance of membrane divided by its capacitance. It is well-known that far from the physiological values of the transmembrane voltage electroporation of cell membranes occurs [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF][START_REF] Michael Leguebe | Conducting and permeable states of cell membrane submitted to high voltage pulses: mathematical and numerical studies validated by the experiments[END_REF], which consists of a high increase of membrane conductance. Moreover the membrane conductance cannot be infinite since it is bounded by the medium conductivity. We thus propose this modified version of the MS model, to account for this phenomenon:

h ∞ (v) = 1 -e -(vgate/v) 2 , 1 τ (v) = 1 τ cl + τ cl -τ op τ cl τ op h ∞ (v), (2.77a
)

∂ t h = g(v, h) = h ∞ -h(v) τ (v) , ( 2.77b
)

I ion (v, h(t, v)) = 1 τ in h(t, v)v 2 (v -1)e -(v/v th ) 2 + v τ out 1 + r max e -(v th /v) 2 , (2.77c)
where g is similar regularisation as to the one given in (2.15), v th v gate is the membrane voltage above which electroporation occurs, and r max 1 stands for the maximal ratio of membrane conductance with membrane capacitance in a fully electroporated membrane. We refer to [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF] for more details.

In the physiological range of membrane voltage, the current behaves similarly to the standard MS model, while for non physiological values of v, the current behaves as a passive conducting pore: I ion ∼ s max v/τ out . The above model is interesting for two reasons. On one side, it provides a simple model for cardiac tissues submitted to electroporation, and on the other side it has a nice mathematical properties and justifies the Lipschitz assumption made earlier on the ionic function, as we will see in the following proposition.

Proposition 17 (Lipschitz property of the modified MS). Let T > 0, and let h 0 be a smooth function such that ∀v ∈ R, h 0 (v) ∈ (0, 1),

|∂ v h 0 | ≤ A.
Then for any v ∈ R, for any t ∈ R + , one has 0 ≤ h(t, v) ≤ 1. Moreover the current I ion of the modified version of Mitchell-Schaeffer model (2.77) is globally Lipschitz with respect to v. That is there exists K such that for any λ 1 , λ 2 , for any t ∈ R + :

|I ion (λ 1 , h(t, λ 1 )) -I ion (λ 2 , h(t, λ 2 ))| ≤ K|λ 1 -λ 2 |, ∀t ∈ R + , ∀(λ 1 λ 2 ) ∈ R 2 .
Proof. By hypothesis on h 0 , as in Lemma 6, the function h stays between 0 and 1 for all t > 0. Then since g is a smooth fonction of v, h is globally Lipschitz with respect to v. Moreover one has the following ODE on H := ∂ v h:

∂H ∂t = 1 τ (v) h ∞ (v) - τ (v) τ (v) (h ∞ (v) -h(v)) -H , ( 2.78
)

H | t=0 = ∂ v h 0 . (2.79) Let à = max v∈R h ∞ (v) - τ (v) τ (v) (h ∞ (v) -h(v)) ,
then H is uniformly bounded on R + by max(A, Ã) and there exists S ∞ such that the function

|∂ v I ion (t, v)| ≤ S ∞ , ∀t ≥ 0. (2.80)
Therefore by definition of I ion one has the global Lipschitz property for any t > 0 and

(v 1 , v 2 ) ∈ R 2 : |I ion (v 2 , h(t, v 2 )) -I ion (v 1 , h(t, v 1 ))| ≤ v 2 v 1 ∂ v I ion (λ, h(t, v 2 ))dλ + h(t,v 2 ) h(t,v 1 ) ∂ h I ion (v 2 , µ)dµ ≤ K|v 2 -v 1 |.

Two-scale convergence

In this section we provide the rigorous justification for the formal derivation given in the Section 2.3. We will follow the idea of two-scale convergence, given in Allaire (1992), [START_REF] Allaire | Homogenization and two-scale convergence[END_REF]. Here we give an overview of a few results we will use in our work.

Definition 18. A sequence of functions u

ε (x) ∈ L 2 (Ω) is said to two-scale converge to a limit u 0 (x, y) ∈ L 2 (Ω × Y ) if , for any function ψ(x, y) ∈ D[Ω; C ∞ # (Y )], we have lim ε→0 Ω u ε (x)ψ(x, x ε ) dx = Ω Y u 0 (x, y)ψ(x, y) dx dy. (2.81)
Theorem 19. From each bounded sequence u ε ∈ L 2 (Ω), we can extract a subsequence that two-scale converges to a limit u 0 (x, y) ∈ L 2 (Ω × Y ). Furthermore, u ε weakly converges in L 2 (Ω) to u(x) = Y u 0 (x, y)dy.

Corollary 20. Let u ε be a bounded sequence in L p (Ω), with 1 < p ≤ ∞. There exists a function u 0 (x, y) ∈ L p (Ω × Y ) such that, up to a subsequence, u ε two scale converges to u 0 .

Proposition 21. Let u ε be a bounded sequence in H 1 (Ω) that converges weakly to a limit u in H 1 (Ω). Then, u ε two-scale converges to u(x), and there exists a function u 1 (x, y) in L 2 (Ω; H 1 # (Y )/R) such that, up to a subsequence, ∇u ε , two-scale converges to ∇ x u(x) + ∇ y u 1 (x, y).

To apply the two-scale convergence to our problem we need to address two problems: first to find the uniform bounds, independent on ε, on functions v ε , u ε , h ε , ∇v ε , and ∇u ε , and secondly to show the convergence of the non-linear part of the problem. In the rest of the chapter we assume that both sub domains, Ω B ε and Ω D ε , are connected.

Hypotheses on the initial conditions

By Corollary 20 and Proposition 21, we infer that

• If for any ε > 0 the initial function v 0 ε (x) satisfies:

v 0 ε L 2 (Ω B ε ) ≤ C, ∇v 0 ε L 2 (Ω B ε ) ≤ C, then v 0 ε (x) 2-scale, L 2 -------→ χ Y B (y)v 0 0 (x), ∇v 0 ε (x) 2-scale, L 2 -------→ χ Y B (y) ∇ x v 0 0 (x) + ∇ y v 0 1 (x, y) .
• If for any ε > 0 the initial value h 0 ε (x) takes values in (0, 1), then we have

h 0 ε (x) 2-scale, L ∞ --------→ χ Y B (y)h 0 0 (x).
From (2.56) taking φ = u 0 ε , and Poincaré's inequality, we have that the initial condition u 0 ε (x) satisfies

u 0 ε L 2 (Ω) ≤ C, ∇u 0 ε L 2 (Ω) ≤ C,
and then

u 0 ε (x) 2-scale, L 2 -------→ u 0 0 (x), ∇u 0 ε (x) 2-scale, L 2 -------→ ∇ x u 0 0 (x) + ∇ y u 0 1 (x, y).
The homogenised problem for (u 0 0 , u 0 1 ) can be then derived, using the test function φ(x) + εφ 1 (x, x/ε), such that φ(x) ∈ D(Ω) and φ 1 ∈ D(Ω, C ∞ # ). We have,

Ω B ε (σ i + σ e )∇u 0 ε (∇ x φ + ∇ y φ 1 + ε∇ x φ 1 )dx + Ω D ε σ d ∇u 0 ε (∇ x φ + ∇ y φ 1 + ε∇ x φ 1 )dx = - Ω B ε σ i ∇v 0 ε (∇ x φ + ∇ y φ 1 + ε∇ x φ 1 )dx.
Passing to the two-scale limit, and after the partial integration, we have the homogenised problem:

-∇ x • Y B (σ i + σ e )(∇ x u 0 0 + ∇ y u 0 1 )dy -∇ x • Y D σ d (∇ x u 0 0 + ∇ y u 0 1 )dy = ∇ x • Y B σ i (∇ x v 0 0 + ∇ y v 0 1 )dy, in Ω,
with the Gauge condition Ω u 0 0 dx = 0, and the correction term:

-∇ y • (1 Y B (σ i + σ e ) + 1 Y D σ d )(∇ x u 0 0 + ∇ y u 0 1 ) = ∇ y • 1 Y B σ i (∇ x v 0 0 + ∇ y v 0 1 ) , in Ω × Y,
with the periodic boundary conditions for y → u 0 1 (x, y) in Y , and the boundary condition:

(σ i + σ e -σ d )(∇ x u 0 0 + ∇ y u 0 1 ) • n Y B = -σ i (∇ x v 0 0 + ∇ y v 0 1 ) • n Y B , on Ω × Γ.
This problem has a unique solution

(u 0 0 , u 0 1 ) ∈ L 2 (Ω) × L 2 (Ω; H 1 (Y )|R).

Uniform bounds

In order to find the uniform bound, we will introduce the energy functional E λ,κ . This is the key idea that is adapted from the work of Allaire [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], and it is suggested in a different context in [START_REF] Collin | Mathematical analysis and 2-scale convergence of an heterogeneous microscopic bidomain model[END_REF]. We thank S. Impériale from M3DISIM Inria team for fruitful discussions. Firstly, we multiply by h ε , v ε , u ε the first, second and third and fourth equations respectively and integrate by parts, thanks to the transmission conditions. We obtain:

1 2 d dt h ε 2 L 2 (Ω B ε ) + Ω B ε g(v ε , h ε )h ε dx = 0, and (2.82
)

1 2 d dt v ε 2 L 2 (Ω B ε ) + Ω B ε I ion (v ε , h ε )v ε dx + √ σ i ∇(v ε + u ε ) 2 L 2 (Ω B ε ) + √ σ∇u ε 2 L 2 (Ω) = 0.
(2.83)

with σ as before,

σ = σ e , in Ω B ε , σ d , in Ω D ε .
Then we infer the following proposition.

Proposition 22. For any λ, κ chosen later on let define the energy functional

E λ,κ (t, v,u, h) := 1 2 e -λt v 2 L 2 (Ω B ε ) + κ h 2 L 2 (Ω B ε ) + t 0 e -λs λ 2 v 2 L 2 (Ω B ε ) + κ h 2 L 2 (Ω B ε ) + √ σ i ∇(v + u) 2 L 2 (Ω B ε ) + √ σ∇u 2 L 2 (Ω) ds.
2. 4 The rigorous mathematical proofs: existence and convergence

Let (v ε , u ε , h ε ) be the solution to problem (2.55). The energy functional E λ,κ satisfies the following equality

E λ,κ (t, v ε , u ε , h ε ) + t 0 e -λs Ω B ε I ion (v ε , h ε )v ε + κg(v ε , h ε )h ε dx ds = 1 2 v 0 ε 2 L 2 (Ω B ε ) + κ 2 h 0 ε 2 L 2 (Ω B ε ) = E λ,κ (0, v 0 ε , u 0 ε , h 0 ε ).
(2.84)

Proof. By definition of E λ,κ (t, v, u, h), using (2.82) (multiplied by κe -λt ) and (2.83) (multiplied by e -λt ) one has

∂ t E λ,κ (t, v ε , u ε , h ε ) + e -λt Ω B ε (I ion (v ε , h ε )v ε + κg(v ε , h ε )h ε ) dx = 0.
Proposition 23. The following inequalities hold:

t 0 e -λs Ω B ε I ion (v ε , h ε )v ε dx ds ≤ C t 0 e -λs v ε 2 L 2 (Ω B ε ) ds, ( 2.85
)

I ion (v ε , h ε ) L 2 (Ω B ε ) ≤ C v ε 2 L 2 (Ω B ε ) . (2.86)
Proof. From the definition of the ionic function 2.14, we see that I ion (0, h(0)) = 0. Then, using the assumption 5 that the ionic function is Lipschitz, one can easily to see that

|I ion (v, h)| ≤ |v|,
0 ≤ h ε (t, x) ≤ 1, |g(v ε , h ε )| ≤ C, ∀t ∈ [0, T ], v ε (t) L 2 (Ω B ε ) ≤ C, I ion (v ε , h ε ) L 2 (×Ω B ε ) ≤ C, ∀t ∈ [0, T ],
and

∇v ε (t) L 2 ([0,T ]×Ω B ε ) ≤ C, u ε (t) L 2 ([0,T ]×Ω) ≤ C, ∇u ε (t) L 2 ([0,T ]×Ω) ≤ C.
Proof. The estimates on h ε (t, x) and g(v ε , h ε ) come directly from the assumption on the initial condition and Lemma 6. Then from the propositions 22 and 23, we have that

E λ,κ (t, v ε , u ε , h ε ) = 1 2 v 0 ε 2 L 2 (Ω B ε ) + κ 2 h 0 ε 2 L 2 (Ω B ε ) - t 0 e -λs Ω B ε (I ion (v ε , h ε )v ε + κg(v ε , h ε )) dx ds ≤ 1 2 v 0 ε 2 L 2 (Ω B ε ) + κ 2 h 0 ε 2 L 2 (Ω B ε ) + t 0 e -λs c 1 v ε 2 L 2 (Ω B ε ) + κc 2 |Ω B ε | ds.
Note that the volume |Ω B ε | < |Ω| < c 3 . Taking λ such that λ 2 > c 1 , and going back to the definition of E λ,κ , we obtain the uniform bound on v ε (t) L 2 (Ω B ε ) , and hence on

I ion (v ε , h ε ) L 2 (Ω B ε )
. Now, going back to (2.83), and integrating over time, we obtain

t 0 √ σ i ∇(v ε + u ε ) 2 L 2 (Ω B ε ) + t 0 √ σ∇u ε 2 L 2 (Ω) ≤ v 0 ε L 2 (Ω B ε ) ≤ C,
which, together with the assumption 4, provides the uniform bounds on the derivatives.

From Ω u ε (t, x)dx = 0 and the Poincare's inequality, we have the bound on u ε L 2 ([0,T ]×Ω) .

Remark 25 (Two-scale limits). Denote with • the extension by zero in the domain Ω D ε . Note that for any function f ε (x) defined on Ω B ε , for which we have a uniform bound

f ε L p (Ω B ε ) < C, it holds that fε L p (Ω) < C. Denote the sequence of solutions (v ε , u ε , h ε )
of the problem (2.55). From the uniform bounds we obtained, by direct application of the theory developed in [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], we have the following convergences

• hε (t, x) 2-scale, L ∞ --------→ χ Y B (y)h 0 (t, x), a.e. t ∈ [0, T ],
• ṽε (t, x)

2-scale, L 2 -------→ χ Y B (y)v 0 (t, x), a.e. t ∈ [0, T ],
• ∇ṽ ε (t, x)

2-scale, L 2 -------→ χ Y B (y)(∇ x v 0 (t, x) + ∇ y v 1 (t, x, y)), • ũε (t, x) 2-scale, L 2 -------→ u 0 (t, x), • ∇ũ ε (t, x) 2-scale, L 2 -------→ ∇ x u 0 (t, x) + ∇ y u 1 (t, x, y), • I ion (v ε , h ε ) 2-scale, L 2 -------→ χ Y B (y)I 0 (t, x), a.e. t ∈ [0, T ], • g(v ε , h ε ) 2-scale, L ∞ --------→ χ Y B (y)g 0 (t, x), a.e. t ∈ [0, T ],
for some functions h 0 (t, •),

I 0 (t, •), g 0 (t, •) ∈ L 2 (Ω), v 0 (t, •) ∈ H 1 (Ω), v 1 (t, •, •) ∈ L 2 (Ω; H 1 (Y B )/R), u 0 ∈ L 2 ([0, T ]; H 1 (Ω)), and u 1 ∈ L 2 ([0, T ] × Ω; H 1 (Y )/R
). These limit functions satisfy the homogenised problem given in formal lemma 1, except the nonlinear functions from the ionic model which are still not identified.

Non-linear convergence

We are now ready to prove that I 0 = I ion (v 0 , h 0 ), and g 0 = g(v 0 , h 0 ). From the given two-scale limits on the functions we have that the equality (2.84) holds for a limit case in the form,

E 0 λ,κ (t, (v 0 , v 1 ), (u 0 , u 1 ), h 0 ) := 1 2 e -λt v 0 2 L 2 (Ω×Y B ) + κ h 0 2 L 2 (Ω×Y B ) + t 0 e -λs λ 2 v 0 2 L 2 (Ω×Y B ) + κ h 0 2 L 2 (Ω×Y B ) + √ σ i (∇(v 0 + u 0 ) + ∇ y (v 1 + u 1 )) 2 L 2 (Ω×Y B ) + √ σ(∇ x u 0 + ∇ y u 1 )
2 L 2 (Ω×Y ) ds.

2.4

The rigorous mathematical proofs: existence and convergence and

E 0 λ,κ (t, (v 0 , v 1 ), (u 0 , u 1 ), h 0 ) + t 0 e -λs Ω Y B I 0 v 0 + κg 0 h 0 dy dx ds = 1 2 v 0 0 2 L 2 (Ω×Y B ) + κ 2 h 0 0 2 L 2 (Ω×Y B ) .
(2.87)

Using the above properties of the non-linear terms, for appropriate choice of λ and κ, we get

E λ,κ (t, v ε -µ ε , u ε -ρ ε , h ε -η ε ) + t 0 e -λs Ω B ε (I ion (v ε , h ε ) -I ion (µ ε , η ε ))(v ε -µ ε ) +κ(g(v ε , h ε ) -g(µ ε , η ε ))(h ε -η ε ) ds ≥ 0.
We develop the above expression, as follows

E λ,κ (t, v ε -µ ε , u ε -ρ ε , h ε -η ε ) + t 0 e -λs Ω B ε I ion (v ε , h ε )v ε -I ion (v ε , h ε )µ ε -I ion (µ ε , η ε )v ε + I ion (µ ε , η ε )µ ε +κ(g(v ε , h ε )h ε -g(v ε , h ε )η ε -g(µ ε , η ε )h ε + g(µ ε , η ε )η ε ) ds ≥ 0.
(2.88)

Furthermore, for the energy functional we have,

E λ,κ (t, v ε -µ ε , u ε -ρ ε , h ε -η ε ) = 1 2 e -λt v ε -µ ε 2 L 2 (Ω B ε ) + κ h ε -η ε 2 L 2 (Ω B ε ) + t 0 e -λs λ 2 v ε -µ ε 2 L 2 (Ω B ε ) + κ h ε -η ε 2 L 2 (Ω B ε ) + √ σ i ∇(v ε -µ ε + u ε -ρ ε ) 2 L 2 (Ω B ε ) + √ σ∇(u ε -ρ ε ) 2 L 2 (Ω) ds.
Hence,

E λ,κ (t,v ε -µ ε , u ε -ρ ε , h ε -η ε ) = E λ,κ (t, v ε , u ε , h ε ) 1 2 e -λt Ω B ε (µ ε (-2v ε + µ ε ) + κη ε (-2h ε + η ε )) dx + t 0 e -λs λ 2 Ω B ε (µ ε (-2v ε + µ ε ) + κη ε (-2h ε + η ε )) dx + Ω B ε σ i ∇(µ ε + ρ ε )(-2∇(v ε + u ε ) + ∇(µ ε + ρ ε )) + Ω σ∇ρ ε (-2∇u ε + ∇ρ ε ) ds.
Now, we use the fact that (v ε , u ε , h ε ) is solution to the problem (2.55), so equality (2.84) holds. Hence, (2.88) becomes,

1 2 v 0 ε 2 L 2 (Ω B ε ) + κ 2 h 0 ε 2 L 2 (Ω B ε ) + 1 2 e -λt Ω B ε (µ ε (-2v ε + µ ε ) + κη ε (-2h ε + η ε )) dx + t 0 e -λs λ 2 Ω B ε (µ ε (-2v ε + µ ε ) + κη ε (-2h ε + η ε )) dx + Ω B ε σ i ∇(µ ε + ρ ε )(-2∇(v ε + u ε ) + ∇(µ ε + ρ ε )) + Ω σ∇ρ ε (-2∇u ε + ∇ρ ε ) ds + t 0 e -λs Ω B ε -I ion (v ε , h ε )µ ε -I ion (µ ε , η ε )v ε + I ion (µ ε , η ε )µ ε +κ(-g(v ε , h ε )η ε -g(µ ε , η ε )h ε + g(µ ε , η ε )η ε ) ds ≥ 0,
Now let us choose, µ ε , ρ ε and η ε as follows

µ ε (x) = v 0 (x) + αφ(x) + εφ 1 (x, x/ε), ρ ε (x) = u 0 (x) + εψ 1 (x, x/ε), η ε (x) = h 0 (x) + αθ(x),
where φ 1 (x, y) and ψ 1 (x, y) are smooth functions that we can choose as close as we need to v 1 (x, y) and u 1 (x, y), respectively. From the construction we have a strong convergences for µ ε , ∇µ ε , ρ ε , ∇ρ ε and η ε . Now, passing to the two-scale limit in the above inequality, we have,

1 2 v 0 0 2 L 2 (Ω×Y B ) + κ 2 h 0 0 2 L 2 (Ω×Y B ) + 1 2 e -λt Ω Y B (αφ) 2 -v 2 0 ) + κ((αθ) 2 -h 2 0 ) dydx + t 0 e -λs λ 2 Ω Y B (αφ) 2 -v 2 0 ) + κ((αθ) 2 -h 2 0 ) dydx + Ω Y B σ i (α∇φ) 2 -(∇ x v 0 + ∇ y v 1 + ∇ x u 0 + ∇ y u 1 ) 2 dydx - Ω Y σ∇(∇ x u 0 + ∇ y u 1 )dydx ds + t 0 e -λs Ω Y B -I 0 (v 0 + αφ) + I ion (v 0 + αφ, h 0 + αθ)αφ +κ(-g 0 (h 0 + αθ) + g(v 0 + αφ, h 0 + αθ)αθ ds ≥ 0,
Finally, we use the equality (2.87), and divide every term by α = 0, to obtain

α Ω 1 2 e -λt φ 2 + κθ 2 + t 0 e -λs λ 2 φ 2 + κθ 2 + σ i (∇φ) 2 dx + t 0 e -λs Ω (I ion (v 0 + αφ, h 0 + αθ) -I 0 )φ + κ(g(v 0 + αφ, h 0 + αθ) -g 0 )θ dx ds ≥ 0,
Then letting α go to zero, we obtain for any functions θ(x), φ(x) t 0 e -λs Ω (I ion (v 0 , h 0 ) -I 0 )φ + κ(g(v 0 , h 0 ) -g 0 )θ dx ds ≥ 0.

2.4

The rigorous mathematical proofs: existence and convergence Thus, we conclude I 0 = I ion (v 0 , h 0 ) and g 0 = g(v 0 , h 0 ).

Main theorem

Theorem 26. Denote with • the extension by zero in the domain Ω D ε . For any time t ∈ [0, T ] denote the sequence of solutions (v ε , u ε , h ε ) of the problem (2.55). Having the two-scale convergences given in remark 25, then (v 0 , u 0 , h 0 , v 1 , u 1 ) is the unique solution of the following two-scale homogenised system:

∂ t h 0 (t, x) + g(v 0 , h 0 ) = 0, in (0, T ) × Ω, |Y B | (∂ t v 0 (t, x) + I ion (v 0 , h 0 )) -∇ x • Y B σ i (∇v 0 (t, x) + ∇ y v 1 (t, x, y))dy -∇ x • Y B σ i (∇u 0 (t, x) + ∇ y u 1 (t, x, y))dy = 0, in (0, T ) × Ω, -∇ x • Y B (σ i + σ e )(∇u 0 (t, x) + ∇ y u 1 (t, x, y))dy -∇ x • Y D σ d (∇u 0 (t, x) + ∇ y u 1 (t, x, y))dy -∇ x • Y B σ i (∇v 0 (t, x) + ∇ y v 1 (t, x, y))dy = 0, in (0, T ) × Ω, with the correction equations -∇ y • σ i (∇v 0 (t, x) + ∇ y v 1 (t, x, y)) + σ i (∇u 0 (t, x) + ∇ y u 1 (t, x, y)) = 0, in (0, T ) × Ω × Y B , -∇ y • [σ e (∇u 0 (t, x) + ∇ y u 1 (t, x, y))] = 0, in (0, T ) × Ω × Y B , -∇ y • σ d (∇u 0 (t, x) + ∇ y u 1 (t, x, y)) = 0, in (0, T ) × Ω × Y D ,
and the boundary and transmition conditions

σ i (∇v 0 (t, x) + ∇ y v 1 (t, x, y)) + σ i (∇u 0 (t, x) + ∇ y u 1 (t, x, y)) • n Y B = 0, on [0, T ] × Ω × Γ, (σ e -σ d )(∇u 0 (t, x) + ∇ y u 1 (t, x, y)) • n Y B = 0, on [0, T ] × Γ, and
h 0 (t, x) = h 0 0 (x), v 0 (t, x) = v 0 0 (x), in Ω, y → u 1 (t, x, y) and y → v 1 (t, x, y) are Y -periodic,

Furthermore, we can recover the classical homogenized and cell equations if we use the relation

v 1 (t, x, y) = N k=1 ∂v 0 (t, x) ∂x k w v k (y), u 1 (x, y) = N k=1 ∂u 0 (t, x) ∂x k w k (y).
Proof. To find the homogenised equations we choose the test function φ(t, x) + εφ 1 (x, x/ε), with φ(t, x) ∈ D((0, T ) × Ω) and φ 1 (x, y) ∈ D(Ω; C ∞ # (Y )). Then by the partial integration and passing to the two-scale limits, using the assumptions on the non-linear parts, we derive the homogenised system of equations.

Remark 27. The homogenised problem given in the main theorem is the parabolic-elliptic version of the homogenised problem we have derived formally in the section 2.3, for

u i ε := v ε + u ε | Ω B ε
, and then u i 0 = v 0 + u 0 , and

u i 1 = v 1 + u 1 .
Remark 28. The existence and uniqueness of the solution of the homogenised (bidomain) model has been proved in [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF]. It is not necessary from conductivity tensors to be symmetric definite positive, although in practice it is assumed.

Discussion

We have proposed a way to model the electrophysiology of the cardiac tissue, that extends the standard bidomain model with the periodic diffusive inclusions. We have a rigorous and practical way to link structural disease in the tissue to macroscopic electrical conductivity in a bidomain model. It can be used in many contexts, including simulation of fibrotic tissue in general and border zones of scars.

From the more practical side, in the following chapters we are going to use high resolution imaging to experimentally obtain insight into the typical geometry/volume fractions in fibrosis.

There are several limitations of the proposed model. The inclusions that we address are purely diffusive, while we can expect to have different kind of cells in these nonexcitable regions. Hence, we neglected the effect we might have from the ionic activity, due to the cells' membrane. Another difficulty that we need to address are transmission conditions given on the interface of the inclusions. We have used standard conditions, usually used in coupling torso and the heart, which might not be appropriate here.

Introduction

We have shown that the micro structure of the cardiac tissue plays a role in AP propagation by affecting the conductivity tensors for the bidomain model. In this chapter we perform a quantitative analysis of these effects. We test different shapes of inclusions in order to build the intuition and learn what kind of changes to expect in the "diseased" tissue. In the next chapter, we perform the simulation on the "real-life" scenario, where we model the AP propagation in the rat heart based on MRI images. We will be able to choose the suitable shape of the inclusions and use the obtained values for the conductivities, from this chapter, without recomputing them. Hence, the results of this chapter can serve as a database for the future experiments on a tissue or organ level. In this way, except for the initial computations that need to be performed on the unit cell, we do not add any computational cost to the simulation of the bidomain model.

For sake of computational simplicity we perform the tests in 2D at first, before going to 3D case. Hence, we solve the cell problems (2.25) -(2.26) and compute the homogenised (effective) conductivities using (2.24). We solve our problems using the standard finite element approach. The variational formulations for the problems (2.25) are given as

Y B (σ i ∇w i j )∇φ i + Γ φ i (σ i e j ) • n = 0, w i j , φ i ∈ W i := w ∈ H 1 (Y B ), w is Y -periodic , ( 3.1)
where i stands for intracellular, and j ∈ {1, 2} i.e. j ∈ {1, 2, 3}, and e j stands for the j-th unit vectors. The variational formulations for the problems (2.26) are given as

Y (σ∇w e j )∇φ e + Γ φ e ((σ e -σ d )e j ) • n = 0, w e j , φ e ∈ W := w ∈ H 1 (Y ), w is Y -periodic . (3.2)
For both, 2D and 3D case we need to set the initial values for the intracellular and extracellular conductivity tensors, σ i and σ e . In the literature there are several proposed values for the conductivities for ventricular tissue. We list some of them in the Table 3.1. While the proposed values are quite different from study to study, we can notice that they are all of the same order of magnitude and that in each case the intracellular and extracellular anisotropy ratios differ significantly one from each other, i.e. the intracellular anisotropy ratio is about σ i l /σ i n ≈ 10, and the extracellular σ e l /σ e n ≈ 2. We perform the simulations with the values given in Rioux (2013), [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF].

Tab. 3.1:

The bidomain conductivities for ventricular tissue, values in [S m -1 ]. The first 5 columns of the table taken from [START_REF] Sundnes | Computing the electrical activity in the heart[END_REF], and the last column from [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF]. Roth (1988) To observe only the effects on the diffusive inclusions on the conductivities we do not vary the directions of the fibers in Y B , and we assume that fibers are aligned with the x-axis. We assume that the conductivity tensor inside of the inclusive diffusion, σ d , isotropic and of the same order as intracellular and extracellular conductivities of the standard bidomain model. Let us note that this is not a necessary assumption from a theoretical point of view, but merely a computational convenience. The non-excitable tissue can be present in the heart in different forms, as explained in the Chapter 1. The electrical properties and structural arrangements in these regions are not well known, hence we have no reason to make an anisotropy assumption. If in the future an anisotropy is observed, our model can easily incorporate this effect as well.

For the numerical simulations we use the software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. For the visualisation we use gnuplot [START_REF] Williams | Gnuplot 4.4: an interactive plotting program[END_REF], scilab [START_REF]Scilab: Free and Open Source software for numerical computation[END_REF] and paraview [START_REF] Ahrens | 36-ParaView: An End-User Tool for Large-Data Visualization[END_REF] software.

Effective conductivities

Settings for 2D cell problems

We deal with the non-dimensional conductivity tensors, with fibers aligned with the x-axis,

σ i = 1.
(σ i ∇w i 1 )∇φ i + Γ φ i σ i 11 n 1 = 0, Y B (σ i ∇w i 2 )∇φ i + Γ φ i σ i 22 n 2 = 0, w i 1 , w i 2 , φ i ∈ W i := w ∈ H 1 (Y B ), w is Y -periodic . (3.4)
The other two on Y ,

Y (σ∇w 1 )∇φ + Γ φ(σ e 11 -σ d 11 )n 1 = 0, Y (σ∇w 2 )∇φ + Γ φ(σ e 22 -σ d 22 )n 2 = 0, w 1 , w 2 , φ ∈ W := w ∈ H 1 (Y ), w is Y -periodic . (3.5)
We compute the new conductivity tensors with the formulas, for j, k = 1, 2,

σ i * kj = σ i kj |Y B | + σ i k1 Y B ∂ y1 w i j dy + σ i k2 Y B ∂ y2 w i j dy , σ e * kj = Y σ kj dy + Y σ k1 ∂ y1 w j dy + Y σ k2 ∂ y2 w j dy , (3.6)
where j ∈ {1, 2} and σ(y) is defined as in (2.12). We perform simulations with simple shapes of diffusive inclusions such as circles, ellipses, and superellipses, ranging the volume fractions of the inclusions from 0.1% to 98%, if possible. More precisely, the unit cell is [-0.5, 0.5] × [-0.5, 0.5] square, with the inclusion generated with the curve

x r x n + y r y n = 1.
(3.7)

Effective conductivities

We use the parametric form of this equation, given with where e = 2 n , and 0 ≤ t ≤ 2π, and the function sgn is a standard sign function,

x =
sgn(t) =        1, t > 0 0, t = 0 -1, t < 0.
(3.10)

We choose to test six types of the inclusions given in the Table 3 As one can notice, in these first six cases we assume the inclusions to be symmetrical with respect to x and y-axes. See Figure 3.1. The difference between the Cases II and III, i.e. between V and VI, is due to the assumption that the fibers are aligned with x-axis. Hence, the elongated inclusions in the Cases II and V are aligned with fibers, and in the Cases III and VI are perpendicular to the fibers.

In the last case, we assume that inclusions are not symmetric with respect to the x-and y-axes. See Figure 3.2. For this case we use e = 1.0, r x = 2r y , and rotate the parametrisation (3.8) for the angle θ = π/4, hence,

x = √ 2 2 (r x cos(θ) -r y sin(θ)), (3.11) y = √ 2
2 (r x cos(θ) + r y sin(θ)).

(3.12)

In this case, Case VII, the volume of the inclusions can be at most 63%, i.e. max |Y D | = 0.63.

We are interested to observe the possible changes in: the speed of propagation, the principal direction of the propagation and in the anisotropy ratios for the intra-and extracellular conductivities.

Case I.

Case IV.

Case II. Case V.

Case III. Case VI. Case VII.

Results in 2D

The results obtained for the cases I to VII are given in the Figures 3.3 The volume fraction for Y D is varied over a range of values depending on the case. In the Table 3.2 we give the maximum possible value of |Y D | for each case.

As mentioned earlier, we are interested in the changes of the anisotropy ratios. In a healthy cardiac tissue, the anisotropy ratios of the intracellular and the extracellular space are

AR int := σ i l * σ i n * , AR ext := σ e l * σ e n * , (3.13)
where σ i,e l * and σ i,e n * are the eigenvalues of the conductivity tensors σ i,e * . For the standard bidomain model, these are conductivities along the fibers and perpendicular to them, i.e. in our case these would be in the directions of x-axis and y-axis. For the modified bidomain model, these directions may change, and we need to diagonalise the conductivity tensors, i.e to find their eigenvalues σ i,e l * and σ i,e n * , and the corresponding eigenvectors as principal directions of the propagation. For all of the first six cases we observe that for the small volume fractions Y D , σ i * and σ e * approach the initial values for σ i and σ e , respectively. It means that when there is no diffusive inclusion, i.e. when there is no inclusions the modified bidomain model is nothing else but the standard bidomain model, with the initial conductivities. On the other hand, for the large volume fractions Y D , we have that σ i * is approaching 0-tensor, and σ e * is approaching σ d . It means that when the inclusion takes up all of the unit cell the tissue is fully non-excitable, and the modified model reduces to the simple diffusive equation with the diffusion tensor of the inclusion. This can be observed in all cases, but it is the most obvious in the Case IV as we are able to test very large volume fractions for Y D . Obviously, the limit is not attainable because the model would not be well posed. Another important observation is that in these six cases we obtain zero values for the off-diagonal terms for new conductivities. As we will see later on, this is due to the symmetry in the shape of inclusions.

While the limit cases and overall behaviour seem to be similar, there are some obvious differences. The dependence on the volume fraction is not linear, and the non-linearity is different for each case. This can also be observed and emphasised when we compute the anisotropy ratios.

In 3.6. We can note that even though the relationship between the old and new conductivity tensors is not a straightforward and obvious, the principal direction of the propagation stays the same. The off-diagonal terms are almost zero, hence we would expect that this shape of inclusions does not have an important effect on the shape of propagation wave, but only on the speed. Similar results are observed in the Case IV. See Figures 3.11 -3.12. Anisotropy ratios do not stay constant and their change resembles the change in the Case I.

Similarities in the way the conductivities and the anisotropy ratios change are observed also among cases II and V, and cases III and VI.

For the Case II and the Case V one can observe that as the volume of the inclusion reaches its maximum, the anisotropy ratios reach their peak. See Figures 3.7-3.8 and Figures 3. 13-3.14. This means that the propagation in the fiber direction becomes much more important than in the perpendicular direction, i.e. the difference becomes much more significant than in the standard model.

For Case III and Case VI we have the opposite situation. See Figures 3.9-3.10 and Figures 3.15-3.16. There is a switch in the principal direction of the propagation that happens when the inclusion is taking more than ≈ 40% , i.e.more than ≈ 50% in the volume. Consequently, we observe a drastic decrease in the intracellular AR from 10 : 1 to 1 : 10.

Finally, the Case VII is somewhat different due to the non-symmetric shape of the inclusion. In this case we observe non-zero non diagonal values for the tensors σ i * and σ e * . To find the principal direction of the propagation in this case we computed the eigenvectors and respective eigenvalues for the tensors σ i * and σ e * , and then we computed the intra-and extracellular anisotropy ratios with the new eigenvalues. The changes in the anisotropy ratios here do not match with any of the previous cases. The results are shown in the Figures 3.17 As we expected, the changes that occur in the conductivity tensors are not simple, and we could not obtain these effects by simply changing the surface to volume ratio in the standard bidomain model. One interesting case that we still did not discuss is the case where we set σ d = σ e . This would mean that the whole inclusion behaves exactly as the extracellular space around the cells in the healthy tissue. In this case σ e * = σ e for all |Y D |, but σ i * would still change depending on the volume fraction and the shape of inclusions. This shows that our problem is indeed a multiscale problem and needs to be treated with care. Black arrows correspond to the eigenvectors for intracellular conductivity tensors, and blue for the extracellular.

Settings for 3D cell problems

The conductivity tensors

σ i =     1.741 0 0 0 0.1934 0 0 0 0.1934     and σ e =     3.904 0 0 0 1.970 0 0 0 1.970     , (3.14)
are scaled like before with ς = 0.1 S m -1 . We assume that the fibers are always aligned with x-axis. We have to solve six cell problems. The first three on Y B ,

Y B (σ i ∇w i 1 ) • ∇φ i + Γ φ i σ i 11 n 1 = 0, Y B (σ i ∇w i 2 ) • ∇φ i + Γ φ i σ i 22 n 2 = 0, Y B (σ i ∇w i 3 ) • ∇φ i + Γ φ i σ i 33 n 3 = 0, w i 1 , w i 2 , w i 3 , φ i ∈ W i := w ∈ H 1 (Y B ), w is Y -periodic . ( 3.15)
The other three on Y ,

Y (σ∇w 1 ) • ∇φ + Γ φ(σ e 11 -σ d 11 )n 1 = 0, Y (σ∇w 2 ) • ∇φ + Γ φ(σ e 22 -σ d 22 )n 2 = 0, Y (σ∇w 3 ) • ∇φ + Γ φ(σ e 33 -σ d 33 )n 3 = 0, w 1 , w 2 , w 3 , φ ∈ W := w ∈ H 1 (Y ), w is Y -periodic . ( 3.16)
We compute the new conductivity tensors with the formulas, for j, k = 1, 2, 3,

σ i * kj = σ i kj |Y B | + σ i k1 Y B ∂ y1 w i j dy + σ i k2 Y B ∂ y2 w i j dy + σ i k3 Y B ∂ y3 w i j dy , σ e * kj = Y σ kj dy + Y σ k1 ∂ y1 w j dy + Y σ k2 ∂ y2 w j dy + Y σ k3 ∂ y3 w j dy , (3.17)
where σ(y) is defined as in (2.12). We perform simulations with the diffusive inclusions in the shape of superellipsoids. More precisely, the unit cell is [-0.5, 0.5] × [-0.5, 0.5] × [-0.5, 0.5] cube, with the inclusion generated with the surface We choose to test four types of the inclusions given in the Table 3.3. In the Case I inclusions take the shape of the smoothed box, and their volume takes up to 98% of the unit cell. In the Case II, we keep the r x = 0.49 fixed, meaning that it almost touches

Effective conductivities

the radius of inclusion is the same in each direction, and for the Case II, where we have fixed r x = 0.49 and vary r y = r z over the range (0.01, 0.49).

In the other two cases, Case III and IV, we assumed that the inclusions are in the shape of plate that varies it's thickness. What we can notice for both of them is that the conductivity in the direction perpendicular to the plate of inclusion drops to the values close to zero even for a small volume of inclusions. This in effect has that the anisotropy ratio σ i 11 * /σ i * 22 , in the Case IV, drops from 9.0 to less than 1.0. In this case we see clearly the importance of the geometry of inclusions.

Discussion

Finally, we can conclude that quantitatively both volume and geometry can play a role in the effective conductivities. Going from 2 to 3 dimensions obviously adds to the complexity of possible shapes of inclusions one could test. We have tested only few possible shapes. The differences between different cases are more obvious when we observe the changes in the anisotropy ratios. The interesting thing is that for a given shapes and volumes of inclusions the principal direction of the propagation changes. Hence, the knowledge of the fiber direction might not be anymore the sufficient parameter to determine the main direction of the AP propagation.

Wave propagation in 2D tissue patch

The next step is to test how these changes in conductivities affect the AP propagation in the tissue. We will perform two kinds of tests: the wave propagation with the simple circular initial condition in the tissue patch, and the spiral wave dynamics. All the tests will be performed on 2D tissue patches.

Settings for homogenised problem

We perform simulations of the standard and modified bidomain model, i.e. we solve the homogenised model given in (2.23). As explained in the Chapter 2, our equations do not depend directly on the small scale, and all of the inclusions effects are included in the conductivity tensors and on the predefined parameter |Y B |. The modified bidomain model reduces to the standard bidomain model in the case |Y B | = 1.0, σ i * = σ i and σ e * = σ e . For the other cases when |Y B | < 1, we directly implement the conductivity tensors we have computed in the previous section, based on the case we want to test. For ionic model we use the standard Mitchell Schaeffer model, see Section 1.2.5. We implement the non-dimensional version of the bidomain model, using the non-dimensionalisation performed in [START_REF] Rioux | A predictive method allowing the use of a single ionic model in numerical cardiac electrophysiology[END_REF], see Section 1.

The non-dimensional version of the model is

∂ t h + g(v, h) = 0, in [0, T ] × Ω, N∇ • σ i * ∇u i = (∂ t v + I ion (v, h))|Y B |, in [0, T ] × Ω, N∇ • (σ e * ∇u) = -(∂ t v + I ion (v, h))|Y B |, in [0, T ] × Ω, ∇u i • n = ∇u • n = 0, on [0, T ] × ∂Ω, v(0, x) = v 0 (x), in Ω, h(0, x) = h 0 (x), in Ω, (3.19)
where the scaling parameters are: V m = 0.125V for potentials u i and u e , and transmembrane voltage, v = u i -u e , that is additionally scaled between 0 and 1; ς = 0.1S/m for conductivities σ i * and σ e * ; T = 10 -3 s for time; L = 5 × 10 -4 m for space; and finally as a consequence we have the non-dimensional parameter N = 0.2. We use the finite element method for the space discretisation with SBDF2 time scheme [START_REF] Ethier | Semi-implicit time-discretization schemes for the bidomain model[END_REF], as in the Algorithm 3. We are solving the equations on the square of a size [-30, 30] × [-30, 30], or 30mm × 30mm. The mesh step is ∆x = 0.16667mm, i.e. 180 × 180 nodes, and the time step ∆t = 0.1ms.

Circular wave propagation

In the first set of tests, we explore the effects of the new conductivities on the propagation of a wavefront in a homogeneous tissue. In other words we assume that the whole patch of tissue has the same properties, i.e. the constant volume fraction of the inclusions and constant conductivities. The initial condition for this set of simulations is given in the Figure 3.28. We assume that the gating variable is h 0 = 0.9 in the whole domain, and the transmembrane voltage is zero everywhere except in the central circle of radius 5mm, where it is v = 0.5.

We run simulations for the standard bidomain model, and plot the snapshots of the transmembrane voltage for the times 20ms, 50ms and 100ms, in the Figure 3.29.

With the same initial condition we perform the simulations for our seven cases given in the Section 3. As we predicted in the cases II and V, (Figures 3.33 and 3.34), we observe emphasised propagation in x-direction or along the fibers, while the propagation in y-direction is faster for the case V than in the case II. This is due to the fact that for |Y B | = 0.61 the inclusion of the shape II almost touches the walls of the periodic cell.

For the cases III and IV, ( Finally for the case VII, (Figure 3.37), we observe a more complex change in the principal direction and the angle of change depends on the volume fraction, as we have already showed in the section 3.2.2.

For all tested cases we can see that the speed of propagation is reduced, although in some cases it is not significant. We can see that when the inclusion is such that it touches the walls of the periodic cell, the principal direction of the propagation moves towards the longitudinal axes of the inclusion.

Spiral wave dynamics

Here we test a more complex AP propagation. The size of the tissue patch that we use it the same as before [-30, 30] × [-30, 30], or 30mm × 30mm, with mesh of 180 × 180 nodes. The time step is dt = 0.1ms. The initial conditions are given as follows: the transmembrane voltage is everywhere zero, except in a small region, -30 < x < -17 and -3 < y < 3, where v 0 = 0.99; the gating variable is zero for y < 0 and h 0 = 0.5 for y > 0. See Figure 3.38.

Firstly, we simulate the standard bidomain model. The results are shown in the Figure 3.39. As one can notice the spiral wave starts to form at the time 80ms, but very soon it disappears as the tissue does not have time to go trough the refractory period.

Let us see what happens when we have the same size of the tissue patch and the same initial conditions, but heterogeneous structure. Part of the tissue is assumed to be healthy where we use the same parameters as for the standard bidomain model. The rest we call the scar and apply different |Y B | and conductivities. We will test three kinds of scars, that we call SCAR L, M and S.

In the first test, we assume that the scar has a circular shape of radius r 2 = 22, in [START_REF] Del Corsso | Transfection of mammalian cells with connexins and measurement of voltage sensitivity of their gap junctions[END_REF][START_REF] Del Corsso | Transfection of mammalian cells with connexins and measurement of voltage sensitivity of their gap junctions[END_REF]. On the left is the scaled transmembrane voltage v 0 , which is zero everywhere except in the small region -30 < x < -17 and -3 < y < 3. On the right is the gating variable h 0 that takes a constant value 0 under x-axis, and the value h 0 = 0.5 above x-axis. The result of this simulation is shown in the Figure 3.42. We can observe that the AP propagation gets blocked by the scar. In the center of the scar the wavefront gets diffused over, and hence the change in the transmembrane voltage v, is not enough to activate the tissue. Hence, the wavefront needs to get around the scar which provides enough time for the tissue that was initially activated to go trough the recovery period and get reactivated once the wavefront reaches it again. We run simulation for the time T = 1.3s and observed the steady spiral around the scar. In the next simulation we use the same shape and radius of the whole scar, but we use a larger border zone, i.e. we reduce the center of the scar to be r 1 = 9 (SCAR M). See Figure 3.40, center. We use the same principle as before of assigning the effective conductivities for our modified bidomain model. See Figure 3.41, top right. The results of this simulation are shown in the Figure 3.43. We run simulation for the final time T = 2s. Like before we observe the steady spiral around the scar. One can notice that the spiral is not going around the whole scar. It means that the volume fraction |Y B | has to be small enough for the AP propagation to be slowed down to obtain the spiral in these settings.

In the third simulation we reduce further the center of the scar to r 1 = 4 (SCAR S). See Figure 3.40, right. The conductivities in this case are shown in the Figure 3.41, bottom left. And the results are shown in the figure 3.44. Here we loose the spiral. The center of the scar is not large enough to slow down the AP propagation in the center of the tissue. We can see though that the shape of the wavefront is changed with respect to the standard bidomain model simulation.

In the last test we go back to the SCAR M, and we change the assumption on the shape of inclusions in the scar. This time we choose the Case VI. See Figure 3.41, bottom right. The results are shown in the Figure 3.45. Interestingly, this time we do not have a spiral. The shape of the wavefront is very distorted. In the center of the scar the AP gets propagated in y-direction faster than in the x-directon, i.e. it gets propagated faster in the direction perpendicular to the assumed fiber direction. This provides the activation of the tissue on the left side of the patch quicker than in the standard case.

Discussion

This section concludes the argument that the shape of non-excitable regions can affect significantly the AP propagation. Of course, our simulations have limitations. We have tested very limited number of shapes for the inclusions, and just a few examples of tissue structure and scar shapes. Already, we can observe the complex changes in the wavefront. Adding more complexity in the shapes of both -inclusions and scars, would give even more complex patterns. The next step in this work can be looking into histological images of the fibrotic hearts, and finding a way to give a better statistical estimation of the tissue heterogeneities. There is a change of the direction of the propagation of AP in the scar, which alternates the shape of the wavefront significantly. This kind of activation does not provide the spiral in the tissue in contrast to the same setting with the inclusions of the type IV. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms, ... , 480ms.

Introduction

In this chapter we present the results of the numerical simulations of the bidomain model and the modified bidomain model on slab of tissue, using the HD-MR Images of the rat heart. The numerical part of this work has been done in collaboration with Yves Bourgault from University of Ottawa, Canada. The HD-MR Images are done by Stephen Gilbert, with whom we also had an important discussions on the topic.

The rat heart ionic model has a much shorter AP duration than a human heart, and the ionic model of AP dynamics is given as a complex physiological model. We use the asymptotic analysis described in the Section 1.3.3, to be able to use Mitchell-Schaeffer model instead and reduce the computational time.

The main part of this chapter is a description of the full process needed to obtain the simulation on the real heart. We start from the HD-MRI of the heart and perform an image analysis to segment the heart, detect fibrotic tissue, estimate its volume fraction, build the mesh for the numerical computations and obtain all the parameters needed to finally perform the 3D simulation.

Modified bidomain model adapted for the rat heart 4.2.1 Equations and dimensions

The problem we solve here is the homogenised bidomain model with modified conductivities, that account for the extracellular space detected on the MR Images. The non-dimensional version is given as follows,

∂ t h + g(v, h) = 0, in [0, T ] × Ω, N∇ • σ i * ∇u i = |Y B |(∂ t v + I(v, h)), in [0, T ] × Ω, N∇ • (σ e * ∇u) = -|Y B |(∂ t v + I(v, h)), in [0, T ] × Ω, ∇u i • n = ∇u • n = 0, on [0, T ] × ∂Ω, v(0, x) = v 0 (x), in Ω, h(0, x) = h 0 (x), in Ω. (4.1)
Here Ω denotes the rat heart. We have to be careful with the parameters scales. To have a correct non-dimensionalisation we need to set the following parameters:

V m = V max -V rest , χ, C m , ς, T and L.
As it is explained in the Section 1.3.3, from here we will have the non-dimensional number,

N = Tς χC m L 2 . (4.2)
First of all let us note that the rat heart has a much smaller size than a human heart. While the human heart is of a size of a fist, measured in centimetres, the rat heart is of a size of the tip of the finger, given in millimetres. Additionally the data that we are working on are given in millimetres, hence we choose L = 0.001m.

The surface to volume ratio, χ, also needs to be adjusted. In the study by Satoh et al, ([96]), they performed the morphological studies on the cardiac myocytes in several species. In the Table 4.1 we show their data obtained in rats. The specific capacitance of cell membranes is approximately C m = 1µF/cm 2 = 0.01F/m 2 , and is fairly constant among different cell types and species ( [START_REF] Hiroshi Satoh | Surface: volume relationship in cardiac myocytes studied with confocal microscopy and membrane capacitance measurements: species-dependence and developmental effects[END_REF]). From here we can recompute the surface to volume ratio as,

χ = C m C m V = 8.44F/l 0.01F/m 2 = 8.44 × 10 5 m -1 = 8.44 × 10 3 cm -1 , ( 4.3)
which is four times larger than in the human heart. The characteristic conductivity, ς, is kept the same as in the human heart. The conductivities of the healthy tissue, σ i and σ e , are kept the same as in the model for the human heart. Based on the literature we have similar conductivities for all mammalian hearts, no significantly different values are given for different species [START_REF] Roth | Electrical conductivity values used with the bidomain model of cardiac tissue[END_REF]. The time scale is free to be fixed and we use T = 0.001s. Now we are able to obtain the non-dimensional number, N = 0.0125.

The rat AP Mitchell-Schaeffer Model

The next step is setting the parameters of the MS model to correspond to the ionic model of the rat heart. The ionic model that corresponds to the AP of the rat heart is given in [START_REF] Sandeep | A mathematical model of action potential heterogeneity in adult rat left ventricular myocytes[END_REF]. As it is explained in the Section 1.3.3, we are not interested in the full expression of the ionic model. We are interested in the duration of each phase of the action potential. The relevant parameters are given in the Table 4.2. We have different values for epicardial and endocardial cells. In order to perform the asymptotic analysis we will use only the values for the epicardial cells.

Tab. 4.2: Action potential characteristics in epicardial and endocardial models, from [START_REF] Sandeep | A mathematical model of action potential heterogeneity in adult rat left ventricular myocytes[END_REF].

Parameter

V We have used the algorithm written by Ngoma, Bourgault and Nkounkou [START_REF] Vianney | Parameter Identification for a Non-differentiable Ionic Model Used in Cardiac Electrophysiology[END_REF], that enables us to automatically determine new parameters of the MS model. We have used the original code, provided by Yves Bourgault, and the details on the algorithm can be found in [START_REF] Vianney | Parameter Identification for a Non-differentiable Ionic Model Used in Cardiac Electrophysiology[END_REF]. Due to the very sharp depolarisation front, and a much shorter AP duration compared to the human heart, some tuning had to be done to avoid the self-excitation in the MS model. The v gate is set to be 0.02, and the time scales for the AP duration and refractory period are adjusted, but the depolarisation time and the total duration of AP signal are preserved to be as in the original model. The resulting parameters for the MS model are The parameters for the bidomain model for the rat heart are summarised in the Finally, the parameters of interest in our study are the modified conductivities, σ i * and σ e * . They are computed from σ i and σ e , given in the Table 4.2, and from the knowledge on the shape and volume of diffusive inclusions, as we have described in the previous two chapters. Also, they depend on the fiber orientation in the tissue. Our goal is to run the simulation of the homogenised model in the real rat heart, based on the HD-MR images, hence we will need to perform some image analysis to precisely define the domain, to obtain the mesh that will be used in numerical simulations and to define modified conductivities.

I ion (v, h) = 1 τ in hv 2 (v -1) + 1 τ out v, g(v, h) =    1 τopen (1 -v), for v < v gate , -1 τ close v, for v ≥ v gate .

Image analysis

In what follows we give a detailed explanation on each step of the image processing we need to perform to be able to run the simulation on the real rat heart. The flowchart on the process of the obtaining parameters from the images is given in the Figure 4.1.

One of the difficulties of this work is dealing with many software and different kinds of files that these software support. In particular, we use the open source software Seg3D [START_REF] Cibc | Seg3D: Volumetric Image Segmentation and Visualization[END_REF], and SCIRun [START_REF]SCIRun: A Scientific Computing Problem Solving Environment, Scientific Computing and Imaging Institute (SCI)[END_REF], developed at the SCI Institute at the University of Utah, to perform most of the analysis. We also use Matlab [START_REF]MATLAB version 8.5.0.197613[END_REF]. Finally, we perform simulations with Freefem++ [START_REF] Hecht | New development in FreeFem++[END_REF], and use Paraview [START_REF] Ahrens | 36-ParaView: An End-User Tool for Large-Data Visualization[END_REF], to visualise our results. As these software do not all support the same type of files, to perform conversions we use few scripts written in the programming language Perl [START_REF] Wall | The Perl programming language[END_REF].

HD-MRI Data

The starting point for us here are the data of the real rat heart. Data provided by IHU-Liryc, Bordeaux1 . The HD-MR Imaging has been performed by Stephen Gilbert, to whom we want to thank for his very insightful comments on the cardiac histology, effects of fiber orientation and fibrotic tissue and detailed explanations of the provided data.

The detailed description of the images can be found in the paper of Gilbert et al 2012, [START_REF] Gilbert | Visualization and quantification of whole rat heart laminar structure using high-spatial resolution contrast-enhanced MRI[END_REF]. The MR Imaging have been performed on the heart of male Wistar rat. After

Segmentation

First of all, we need to have a defined domain of simulation. For this we use the segmentation techniques to extract the region of interest from the MR Images. We use Seg3D software to mark the slab of tissue we are interested in. We use several filters to mark the region of interest. When performing segmentation many things have to be corrected manually. One should be careful that the region of interest is connected, that the filters did not produce an accidental holes in the tissue. Additionally one should be careful with the boundaries and make sure that there are not too small or too sharp corners, which could affect the mesh generator to produce "bad" triangles in the mesh. We have used several filters and tools available in the Seg3D software:

• Median filter -it is a non-linear filter that reduces noise and preserves or enhances sharp boundaries. For every pixel (except the edges), this filter finds the median value of the grey scale of the neighbouring pixels and uses this value to replace the original pixel. The neighbourhood is controlled by the distance parameter. We use the distance parameter that equals 3.

• Threshold -the tool that is always applied on the data layer, and creates the mask.

The mask is the image of the same size where each voxel has assigned values 0 or 1, based on some criteria. When we use the tool threshold, it is based on the range of the grey scale values that we want to segment. In other words, this tool is used to split roughly voxels into those who belong to the heart tissue (that get value 1) and those that belong to the void (that get value 0).

• Paint brush -the tool used to fix the noise. After the threshold is applied on a data layer, there are some voxels that get assigned wrong values. This creates unrealistic holes in the heart tissue, as well as some artificial "heart" voxels around the heart borders or in the void. Hence, we use the paint brush manually, layer by layer, to fix "bad" voxels.

• Crop -the tool for assigning values 0 or 1 based on the position of the voxel in the image. We use this tool to determine the part of the heart that we want to run our simulation on. We crop out the part of the wall of the left ventricle.

• Arithmetic filter -runs basic arithmetic and boolean operations on both image and mask data on a voxel by voxel basis. We use the boolean AN D to combine the mask crated by the tool crop and the mask that represents the heart tissue, to obtain the correct domain for our simulation.

The result of our segmentation is shown on the Figure 4.3. One should notice that even though the Seg3D software works with the image layers and voxels, the output files preserve the information on the space scaling, so we are able to pass them to some other software and work with metric dimensions.

Defining diffusive inclusions

The images that we use have been done by Stephen Gilbert and collaborators in order to study the laminar structure of the heart tissue. In our case we use these images to simulate the AP propagation in the fibrotic tissue. There is an obvious disadvantage, as the presented heart is healthy in a real-life. The gaps between layers are enlarged with the perfusion for the imaging purposes, to emphasize the myofiber and myolaminar structure in the healthy heart. In reality these gaps are much smaller. With this in mind, we are making a modelling assumption: that these gaps represent the fibrotic or fatty tissue that we interpret as diffusive inclusions. Our goal is to demonstrate that based on HD-MRI we are able to reconstruct and implement the modified bidomain model with a reasonable approximation. Once we have the relevant MRI we would be able to obtain more accurate modelling of the inclusions in a heart.

The next step is to decide what is considered to be the diffusive inclusion. We need to go back to the original HD-MRI to get into the smaller scale image analysis. On the MRI the tissue is represented with white and the extracellular space with black voxels. The simplest way to decide if the voxel belongs to the "normal" cardiac tissue or to the inclusion is by using the threshold on grey scale in MRI. Since we work on a small scale, we do not perform any additional fixing of the shapes of the inclusions, as it would be daunting work and we have no information based on which we would improve the results manually. Hence, we are satisfied with a simple threshold approach. After this we use the Arithmetic filter to extract the inclusions only on the segmentation, and not on the complete image. See

Image analysis

We compare different thresholds in one layer of the segmentation in the Figure 4.5. And the closer look for the first threshold value is shown in the Figure 4.6.

Computing volume ratio: |Y

B | = 1 -|Y D |
In the modified bidomain model we proposed there is an effect of diffusive inclusions on the conductivities, while both volume ratio and the geometry of inclusion play a role. Here on the other hand, we assume that inclusions are all of the same shape, the shape of super-ellipsoid whose surface is given with the equation

x r x n + y r y n + z r z n = 1.
This approach allow us to be able to use the results from the Section 3.2.4, where we take the results for the case I, i.e. when r x = r y = r z . In reality the inclusions are not all of the same shape, and we should perform some additional image analysis to assess the best approximate shapes. While this would be a better approximation, we keep things simple and use only this type of inclusions throughout the tissue. Now if we have the geometry of the inclusions, we focus on analysing their local volume ratios. For this purpose we use Matlab software. The input file is the matlab file (exported from Seg3D), that contains image of the size 256 × 256 × 512, where voxels have value 1 if they belong to the inclusion on the segmentation, and all others are zero. We apply the following algorithm:

1. Around each voxel of the image use the window of a given size, we work with the window size W := [5 × 5 × 5]; 2. Count the number, num, of the voxels inside of the window that belong to inclusions, here 0 < num < 125. 3. Devide num with the total number of voxels in the window, here with 125. As one can notice, the local volume ratio will depend on two parameters: 1) the threshold for the inclusion detection, and 2) the window size. What we perform with Matlab is similar to the Seg3D mean data filter. We could not use Seg3D in this case, because we had to work with the values on the voxels in the mask, and not on the data layer, which in our version of the software is not supported. For the full list of the available filters in Seg3D and the detailed explanations we refer the reader to the wikipage: http://scirundocwiki.sci.utah.edu/SCIRunDocs/index. php5/CIBC:Seg3D2:Reference .

In the case of lower pixel density of the image, for example in the case of bigger hearts or for the lower resolution of the MRI, one could directly use the relative value of the grey scale intensity to decide the volume fraction. Say, maximum value of the grey scale in the image is M , and minimal is m. If the grey scale of the given voxel is x, the volume fraction for this voxel would be:

x-m M -m .
image, onto nodes of the mesh. Let us note that both files contain the information on the physical size and coordinates of the voxels, that is used for the correct mapping of the information.

The volume fraction mapped on the mesh is shown in the Figure 4.8.

Fibers directions on mesh

As we have described in the introduction, section 1.3.2, the fiber structure in the cardiac tissue leads to anisotropy in the bidomain model. Hence, for the simulation of the bidomain model it is important to assess the fiber orientation in the tissue.

Various techniques can be used to obtain the fiber direction such as histology, optical techniques and diffusion tensor MRI (DT-MRI). One way to impose the fiber directions in the ventricles is rule-based fiber orientation, where the direction changes smoothly from endocardium to epicardium [START_REF] Knisley | Roles of electric field and fiber structure in cardiac electric stimulation[END_REF]. This is a reasonable approximation in a case where we do not have any access to the realistic fiber orientation.

The more precise technique is use of the DT-MRI. It is a non-invasive imaging method that uses the diffusion of molecules to generate a contrast on MR images. This technique exploits the fact that the diffusion of molecules is not free, but it is constrained with obstacles such as cell membranes. For our study, we will use the data from DT-MRI.

The provided data from DT-MRI are given on the four times coarser image than the original HD-MRI data, i.e. 64 × 64 × 128, but for the same view field, i.e. the spacing between voxels is four times bigger, and it is 0.2 × 0.2 × 0.2mm. The file is given in .VTK format, and it contains vectors of the fiber orientations. To be able to use these data in our simulation, we have to map it on the mesh nodes. For this purpose, we use again SCIRun software. Our version of the software does not support the import of the vector files in .VTK format, so we had to split the given file into three files that contain separately information for x-, y-and z-coordinates of the vectors. Then we use the scheme in SCIRun to combine them back into one file of the .NRRD format2 , and to visualise them as in the Figure 4.9, left. As one can notice, we have fiber orientation provided on the full heart. To map fibers on the segmentation and on the mesh, like before, we use MapFieldDataOntoNodes. The results are presented in the Figures 4.9, right, and 4.10. The software gives the continuous bidirectional vector field throughout the domain. 2012 [START_REF] Gilbert | Visualization and quantification of whole rat heart laminar structure using high-spatial resolution contrast-enhanced MRI[END_REF], the original data on the full heart (left) and the projection of the fibers on the segmentation (right). Visualised with SCIRun.

Fig. 4.10:

The bidirectional vector field of the fiber orientations projected on the mesh nodes, visualised with SCIRun. This field is much finer than the one in the Figure 4.9, due to the much higher number of nodes in the mesh, compared to the number of voxels in DT-MR image.

Simulation in FreeFem++

At this moment we have obtained all required data from the images. We have segmented and meshed the domain of interest, obtained the volume fraction of the inclusions throughout the domain, and finally we mapped the fiber orientations and volume fraction on the mesh nodes. All these informations are saved in the .VTK files, that are the input files to our FreeFem++ code for running the simulation. Now we focus on the numerical scheme and the additional processing that needs to be done in the FreeFem++ program.

Effective conductivities based on the HD-MRI

Before proceeding to the final stage of the simulation, where we can apply the numerical scheme and solve the system of PDEs, we need to make sure that the conductivities are properly set on each node of the mesh. Perl is used to convert .VTK3 files into .MSH4 file for the mesh, and simple text files for the volume fractions and fiber directions. This part of the work is performed in Freefem++.

The intracellular and extracellular conductivity tensors for the standard bidomain model are constant throughout the heart tissue and their scaled eigenvalues, σi l,n and σe l,n , are given in the table 4.3. We made an assumption that all diffusive inclusions are of the same shape, a shape of superellipsoid as in the Case I, described in the section 3.2.3. In the section 3. Note that these conductivities still do not include the information on the fiber directions, i.e. they are given as if all the fibers were aligned with the x-axis. So, additional processing needs to be done to obtain the conductivities that we will use directly in the simulation. For each node of the mesh we have stored the direction of the fiber, as a vector a l = (f x , f y , f z ). Then, based on section 1.3.2 and equation (1.63), we compute the final conductivities as

σ i * =     σ i * n + (σ i * l -σ i * n )f x f x (σ i * l -σ i * n )f x f y (σ i * l -σ i * n )f x f z (σ i * l -σ i * n )f x f y σ i * n + (σ i * l -σ i * n )f y f y (σ i * l -σ i * n )f y f z (σ i * l -σ i * n )f x f z (σ i * l -σ i * n )f y f z σ i * n + (σ i * l -σ i * n )f z f z     ,
and the same for the extracellular conductivity. Note that for the standard bidomain model, when computing the conductivities, the starting point are the values given in the table 4.3. And we only need to apply the effect of the fibers, since the volume fraction is everywhere

|Y B | = 1.

Space and time discretisation

Now we have all the required parameters for the computation. To run the simulation we need a time-space numerical scheme for solving the modified bidomain model 4.1. The modifications are accounted in the conductivities, σ i * and σ e * , and the parameter |Y B |, otherwise the homogenised model is exactly the standard bidomain model. The proof for existence and uniqueness for the bidomain model, for the isolated heart (without torso), can be found in [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF][START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF][START_REF] Colli | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF]. We are free to proceed directly to the numerical solutions of our problem.

For the space discretisation we use the finite element approach. We need the appropri-ate variational formulation of the problem. We use the one of the parabolic-elliptic problem, as follows: for all t > 0, find the solution (v(t, •), u e (t, •), h(t, •)) ∈ H 1 (Ω)×H 1 (Ω)×L ∞ (Ω), such that Ω u e = 0, and

|Y B | Ω (∂ t v + I ion (v, h)) φ + N Ω σ i * ∇v ∇φ + Ω σ i * ∇u e ∇φ = 0, Ω σ i * + σ e * ∇u e • ∇ψ + Ω σ i * ∇v • ∇ψ = 0, ∂ t h + g(v, h) = 0, (4.6)
where φ(t), ψ(t) ∈ H 1 (Ω) and Ω ψ(t) = 0. The initialisation of the AP propagation is achieved by setting the proper initial conditions v 0 and h 0 . For the time discretisation for the bidomain model we use the second order semiexplicit backward differential formulae (SBDF2), [START_REF] Ethier | Semi-implicit time-discretization schemes for the bidomain model[END_REF]. Let N be a given integer and take the uniform partition of the interval [0, T ] into intervals of the size ∆t = T /N . We look for the numerical solutions of our problem at the times t n = nT /N = n∆t, where n ∈ 0, 1, .., N . Denote by (v n , u e n , h n ) the approximation of (v, u e , h) obtained at the time t n . For the first step we need to apply the Explicit Euler scheme, to obtain the solution at the time t 1 , and to be able to apply SBDF2 for t n , n ≥ 2. See the algorithm 4.

Results

For our simulations we use the FreeFem++ software. We run one test case with several implemented solvers such as conjugate gradient (CG), MUMPS, sparse solver, etc. Finally we use CG solver as it has the most suitable performance from the memory and speed point of view.

We perform the first simulation on the standard bidomain model, including the information on fiber orientation and with the following parameters:

• time step, ∆t = 0.05,

• final time, T = 500,

• initial data: v 0 = 0.9e -2((x-9.5) 2 +(y-8) 2 +(z-6) 2 ) , and h 0 = 0.9.

We use Paraview( [START_REF] Ahrens | 36-ParaView: An End-User Tool for Large-Data Visualization[END_REF]) to visualise the results of the simulation, As we have expected we have the reduction in the velocity of the propagation, but not on the shape of the wavefront. This change in the velocity does not seem to be significant. One should keep in mind that data from HR MRI are on the healthy heart, hence we could not observe significant fibrosis in the tissue. Note also that we used the values of σ d that are of the same scale as the intracellular and extracellular conductivities. Hence, one might expect different results if the value of σ d is much smaller than that. Algorithm 4 Discretisation scheme for modified bidomain model 1: Define ∆t, T and N := T ∆t . 2: Define var. form. space H 1 (Ω) with P1 elements. 3: Declare u i n , u n , v n , h n , φ, ψ ∈ H 1 (Ω), for n = 0, .., N. 4: Define v 0 := v 0 , h 0 := h 0 . 5: Solve the coupled system with EE to obtain v 1 , u 1 :

|Y B | Ω v 1 φ -|Y B | Ω (v 0 φ + ∆tI ion (v 0 , h 0 )φ) +∆tN Ω σ i * ∇v 1 ∇φ + Ω (σ * ∇u 1 ) ∇φ + Ω σ i * ∇v 1 ∇ψ + Ω (σ i * + σ * )∇u 1 ∇ψ = 0 6: Solve the ODE to obtain h 1 : Ω h 1 φ = Ω (h 0 + ∆tg(v 0 , v 0 )) φ. 7: for n = 1 to N do 8:
Solve the coupled system with SBDF2 to obtain v n+1 , u n+1 :

|Y B | Ω 3 2 v n+1 φ -|Y B | Ω (2v n - 1 2 v n-1 ) + ∆t (2I ion (v n , h n ) -I ion (v n-1 , h n-1 )) φ +∆tN Ω σ i * ∇v 1 ∇φ + Ω (σ * ∇u 1 ) ∇φ + Ω σ i * ∇v 1 ∇ψ + Ω (σ i * + σ * )∇u 1 ∇ψ = 0 9:
Solve the ODE to obtain h n+1 :

Ω 3 2 h n+1 φ = Ω 2h n - 1 2 h n-1 + ∆t (2g(v n , h n ) -g(v n-1 , h n-1 )) φ. 10:
Define v n+1 := u i n+1 -u n+1 .

Discussion

What we demonstrated in this chapter is that our model is applicable in real-life settings, i.e. for the heart for which enough imaging data are available. We have used the HR MRI data for the rat heart. Since this heart is very small the images contained a detailed structure of the ventricular myocardium, so we were able to define the fibrotic inclusions. We have segmented a part of the left ventricle of the heart, and analysed images to obtain the local volume fractions of the fibrotic tissue. With this, we run the simulation of the modified bidomain model, without additional computational cost compared to the standard bidomain model. We conclude that we have a rigorous and practical way to link the effects of the structural changes in the tissue, as observed in HR MRI, to the changes in the macroscopic electrical behaviour. In the future work we would like to use our model to asses the modified conductivities in diseased tissues, especially in such settings where the effects would be more visible (e.g in scars). difference of membrane potential between adjacent cells.

We propose a 2D/3D microscopic mathematical model that takes into account the non-linear behaviour of GJCs. We use non-standard geometrical setting for the tissue micro structure, where we abandon the physical connection between adjacent cells. We keep the functional relation between the intracellular potentials between the adjacent cells via boundary conditions.

There has been some interesting work done on the microscopic modelling of the tissue, on a scale of several tenths to a couple of hundreds cells [START_REF] Stinstra | Incorporating histology into a 3D microscopic computer model of myocardium to study propagation at a cellular level[END_REF][START_REF] Jeroen | On the passive cardiac conductivity[END_REF][START_REF] Jeroen | Modeling the passive cardiac electrical conductivity during ischemia[END_REF]. The model proposed in these works is very interesting, as their simulation is given in 3D, assumes a complex shapes of myocytes which are more realistic than a simple rectangles, and provide a good insight on the effects of the parameters. In their study they show that different conductance of the gap junctions leads to quite different electrical behaviour of the tissue, including the effect on the average conductivity and anisotropy ratio. On the other hand, they model gap junctions as a simple resistors, and do not incorporate dynamical behaviour of the GJCs that we are interested in.

Experimental Data on Cx43 and Cx45 in rat

The technique used to perform the dual voltage clamp experiment is described in details in del Corso et al. (2006). This technique was first proposed by [START_REF] Harris | Kinetic properties of a voltage-dependent junctional conductance[END_REF]. It is used to assess the behaviour of GJCs between a single pair of adjacent cells. The patch clamp electrodes are inserted in each cell and a constant transjunctional voltage V j is applied, which provides the insight into the time evolution of the transjunctional current I j , and junctional conductance G j = V j /I j . See Section 1.4.

The behaviour of GJCs is described with the following parameters: G j,max and G j,min represent the maximal and minimal junctional conductances, respectively, and G j,0 and G j,∞ represent instantaneous and steady state junctional conductances. Due to limitations of the experimental setup, G j,max and G j,0 are not necessarily the same. The dynamics is described by τ j,∞ that describes the time rate at which the GJCs transit from the instantaneous to the steady state. The conductivities G j,0 and G j,∞ , and τ j,∞ are measured with respect to transjunctional voltage, V j . We use G j,0 , as a reference value for the maximal conductance that can be achieved in alive tissue. See the section below for details.

All parameters depend on the type of connexins that form GJCs. In healthy ventricular myocytes the most expressed connexin is Cx43, and in much smaller amount Cx45. Hence we focus on GJCs formed by these two connexins. The experimental data (Desplantez et al. 2004[START_REF] Desplantez | Gap junction channels and cardiac impulse propagation[END_REF][START_REF] Desplantez | Influence of v5/6-His tag on the properties of gap junction channels composed of connexin43, connexin40 or connexin45[END_REF], shown in Figure 5.1, suggest that the conductance of homomeric homotypic Cx43 and Cx45 GJCs is symmetric with respect to V j = 0 mV. While we know that GJCs can express significantly different behaviour in case of heteromeric or heterotypic GJCs, for the modelling purposes we will focus only on these two types of the GJCs.

Non-linear gating model

Even though the behaviour of GJCs has been for a long time observed and described as non-linear (Harris 1981), currently the common way to model GJCs in the tissue is to assume passive behaviour. This means that the electrical current for a given transjunctional voltage V j ,

I j = G j V j , ( 5.1)
is assumed to be linear, i.e the junctional conductance G j is assumed to be constant. This simplification is useful and in the healthy tissue it is a good approximation, due to the very different time scales for GJCs (order 3s) compared to ionic channels (order 20-300ms).

Here we want to observe the importance of the dynamical behaviour of GJCs. Hence, we use non-linear functions to fit the experimental data presented above. We assume that GJCs behave as gates and the data are used to specify the nature and behaviour of the gating variables. The rates at which the gates open or close depend on the transjunctional voltage. We write the conductance as follows

G j (t, V j ) = G j,max g j (t, V j ), (5.2)
where G j,max is the above mentioned maximal junctional coupling, and g j represents the non-linear gating variable, that takes values between 0 and 1. We use the work of [START_REF] Harris | Kinetic properties of a voltage-dependent junctional conductance[END_REF] [START_REF] Spray | Equilibrium properties of a voltage-dependent junctional conductance[END_REF][START_REF] Harris | Kinetic properties of a voltage-dependent junctional conductance[END_REF], to describe the model for the gating variable. The reasoning is similar to this described in Section 1.2.3. The channels can be in one of two states: open or close,

C α β O, ( 5.3)
and α and β are rate constants that are voltage dependent. The difference with the standard ionic channels, is that we assume that the closed state is not completely closed, and we can refer also to these two states as high-and low-conductance states. When the voltage is kept on the constant value, α and β are also constant. Then, if n is the fraction of open channels, the dynamics of n at the given constant voltage is given as

dn dt = -βn + α(1 -n). (5.4)
The solution is,

n = n ∞ + (n 0 -n ∞ )e -t τ j,∞ , (5.5)
where n 0 and n ∞ are the initial and the steady state values of n, respectively, and τ j,∞ is the time constant, such that

n ∞ = α α + β
, and, τ j,∞ = 1 α + β .

(5.6)

Non-linear gating model

As we will see further, the parameters g j,∞ and τ j,∞ , can be fitted from the experimental data. Finally, using (5.7) the gating variable dynamics is expressed as

dg j dt = g j,∞ (V j ) -g j τ j,∞ (V j ) , ( 5.16)
with the solution for a constant V j , g j = g j,∞ + (g j,0 -g j,∞ )e -t τ j,∞ .

(5.17)

Data fitting

The dynamics of the gating variable is described by a simple ODE (5.16), where g j,∞ and τ j,∞ are fitted to the experimental data (Figures 5.1 We fit the relationship between g j,∞ and V j to a Bolzmann equation 5.14 where V j,0 is the voltage at which the conductance is half-minimal, g j,min is a normalised voltageinsensitive residual conductance and A is a parameter defining the steepness of voltage sensitivity. We fit the relationship between τ j,∞ and V j in an exponential way, as in (5.15). The parameters obtained from the data fitting in the case of homotypic Cx43 and Cx45 GJCs are presented in Table 5.1.

0D simulation

We perform 0D simulation for the gap junction current with the non-linear model for the GJ conductance, i.e. I j = G j,max g j (t, V j )V j , (5.18)

g j (t, V j ) = g j,∞ (V j ) + (g j,0 -g j,∞ (V j ))e - t τ j,∞ (V j ) , (5.19) g j,∞ (V j ) = 1 -g j,min 1 + e A(V j -V j,0 ) + g j,min , (5.20) τ j,∞ (V j ) = ae -bV j , ( 5.21)
where the parameters G j,max , g j,min , A, V j,0 , a, and b, are obtained from the data fitting, and are given in the Table 5.1. Initially all the channels are assumed to be open, hence g j,0 = 1. Note that we do not use the fit function for the instantaneous conductance.

In the Figure 5.3 we show the results obtained in the study of [START_REF] Desplantez | Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties[END_REF], from which we obtained the data for the fitting of the g j,∞ and τ j,∞ . In the Figure 5.4, we plot the currents for fixed V j = ±10, ±40, ±70, ±100, and ±130 mV. The numerical results for Cx45 correspond to the experimental observations, i.e. are of the same order. As for the results for Cx43 we can observe that the model gives the currents that are almost ten times higher than those in the experiment. This can be explained in the fact that G j,max = 68.3 nS. This value is obtained from another experiment, and as we already discussed in the introduction (Section 1.4.2), the limitation of the patch clamp experiment leads to underestimation of the gap junction current when the maximal conductance is larger than ≈ 10 nS.

In order to recover the experimental results, we reduce the maximal conductance The main reference in this section is the work of Beauchamp et a.l 2012, [START_REF] Beauchamp | Electrical coupling and propagation in engineered ventricular myocardium with heterogeneous expression of connexin43[END_REF]. In this study they quantify the relation between the degree of heterogeneity in Cx43 expression and disturbances in electric propagation. The full details of the experiment can be found in [START_REF] Beauchamp | Electrical coupling and propagation in engineered ventricular myocardium with heterogeneous expression of connexin43[END_REF]. They engineered patterns of murine strands (4-5 mm in length and 50, 100 or 200 µm in width) from various mixtures of wild-type Cx43 that co-express Cx43 and low levels of Cx45, and Cx43KO cells (knockout of Cx43, i.e. cells that do not express Cx43).

The propagation velocity was measured from the time difference in average activation and the distance between 2 regions of interest in the various cell mixtures. The original velocity of 30.5 cm/s in the strands with only wild-type Cx43 cells, dropped to 76%, 55% and 19% of the original velocity when the ratios of Cx43 and Cx43KO cells were 80%:20%, 50%:50% and 20%:80% respectively. Finally the velocity drops to 2.1 cm/s for the strands with only Cx43KO cells.

We would like to develop a model that replicates such behaviour. We would like to study if this change in dynamics could be seen as the direct consequence of the change in the gap junction type. For the modelling purposes we assume that the Cx43KO cells behave as the homotypic Cx45 channels.

Model for a strand of cells

We propose a microscopic mathematical model for the signal propagation in the tissue. We use a non-standard geometrical setting for the tissue micro structure. Namely, we represent each myocyte as a separate rectangular domain, where the boundary represents the cell membrane, see Figure 5.6. Furthermore, we split the boundary into "ionic" and "junctional" boundaries and we apply the ionic model on the former and the previous nonlinear GJ model on the latter. In this way, even if we do not have the physical connection, as is the case in the standard settings, we keep the functional relation between adjacent cells via boundary conditions. The governing equations for the potentials in both intra and extracellular spaces are Laplace equations.

Ω e Γ gr Γ app Ω i Ω i Γ j Γ j
Γ ion Ω e Fig. 5.6: Geometric settings for 2D mathematical model. Ω i , Ω e -intra and extracellular domains, Γ ion -"ionic" boundary, Γ j -"junctional" boundary.

The ionic properties of the membrane are taken into account as the boundary conditions on the "ionic boundary".

∂ t V m + I ion (V m , h) = -σ i ∇u i • n, ∂ t V m + I ion (V m , h) = -σ e ∇u e • n, ∂ t h = f ion (V m , h),        on Γ ion .
Here, V m = u i -u e is transmembrane potential, σ i and σ e are intra and extracellular conductivities and h is a vector of variables in the ionic model. For this study we choose Beeler Reuter ionic model, [START_REF] Go | Reconstruction of the action potential of ventricular myocardial fibres[END_REF]. The non-linear gap junctions are taken into account as a boundary condition on the "junctional boundary". G j,max g j V j = -σ i ∇u i • n, ∂ t g j = (g j,∞ (V j ) -g j )/τ j,∞ (V j ). on Γ j .

(5.22)

Here, V j represents the transjunctional jump of potential in the adjacent cells, g j,∞ is given as (5.14) and τ j,∞ as (5.15). We assume that the system is isolated, i.e. we have the homogeneous Neumann boundary conditions on the external boundary, and on Γ j , σ e ∇u e • n = 0, on ∂Ω e \ Γ ion and Γ j .

(5.23)

We set the initial conditions such that the system is in the steady state. In order to observe the signal propagation we need to provide a stimulus. This is done as follows: on one side we keep the "ground" and on the other we apply the difference of potential, V stim , that can trigger AP in the first cells that are affected. In other words, during the stimulation time t ∈ [t 0 , t 0 + t stim ], where t 0 is the beginning of stimulation, we solve the system with Dirichlet boundary conditions u e = 0, on Γ gr , (5.24) u e = V app , on Γ app , (5.25) G max g j (V j )V j φ i = 0.

Here we do not deal with the well-posedness of the problem. This is in fact the modified microscopic problem from [START_REF] Colli | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF][START_REF] Veneroni | Reaction-diffusion systems for the microscopic cellular model of the cardiac electric field[END_REF]. In their work they prove the existence and uniqueness of the microscopic problem for some membrane ionic model. The difference is in the setting of the geometry, as they assume intracellular space to be be connected. And they do not consider gap junctions as the boundary condition, but rather as physical pathways between adjacent cells. In such settings it is not obvious that the non-linear dynamics of gap junctions can be considered.

Preliminary numerical analysis

To perform the numerical simulations we use the finite element approach with semiimplicit time discretisation scheme. We split intracellular and extracellular problems and use the iterative scheme for each time step.

In this 2D setting, for each cell we discriminate between left and right gap-junctional boundary that we denote Γ j,l and Γ j,r , respectively. We have defined the transjunctional potential as V j (x) := u i (x -G) -u i (x) on Γ j,l , or V j (x) := u i (x) -u i (x + G) on Γ j,r . Then we split the intracellular and extracellular problems, to solve with the iterative scheme. The two problems are the following: (u e +V m -∆tI ion (V m ), h)v i -Γ j,l ∆tI j v i + Γ j,r ∆tI j v i = 0.

We have tried this approach with 4 and 10 myocytes, for which the CPU time to simulate 10 -20ms is of order of 10hrs. The result on 4 myocytes is shown in the Figure 5.7. We have also tried to change the boundary conditions on the outside domain to be fully Neumann, and to apply external stimulus as the source term in the extracellular space. The results were similar. At the moment we did not make any conclusion on the observations. The next thing we did is to investigate the behaviour of our model by assuming that there are no gap junctions, i.e. we would apply the ionic model on both Γ ion and Γ j . The iteratvie semi-implicit scheme we used is described in the algorithm 5. The size of the myocytes used in simulations are 100µm × 20µm. This correspond to the average length and width of human ventricular myocytes. The space between adjacent cells is set to 1µm. The mesh step is 1µm and the time step in the simulation is 0.02ms.

In the first numerical experiments we used only one myocyte embedded in the extracellular space. By stimulating extracellular space we were able to trigger the AP in the cell and observe its propagation as expected. We observed all of the four phases of the Algorithm 5 Iterative scheme for microscopic model without GJ 1: Define ∆t, T . 2: Define threshold. 3: Define space V hi on Ω i with P1 elements, and declare u i , u 1 i , v i ∈ V hi. 4: Define space V he on Ω e with P1 elements, and Dirichlet boundary conditions on Γ gr and Γ app , and declare u i , u 1 i , v i ∈ V hi 5: Define space V h on Ω with P1 elements, and declare V m and all the state variables of BR model, h, to be in V h. They take 0 values everywhere except on the boundary Γ ion . 6: Define V m , h on Γ ion , and V j on Γ j . At this moment u i (t) and u e (t) converged to certain values and we compute: V m (t) = u i (t) -u e (t) on Γ ion

22:

Define all the variables of BR model on Γ ion accordingly.

23:

Compute: I ion = I ion (V m (t), h).

AP in a single cell simulation. On the other hand in two-cell experiments we were not able to observe the propagation of AP from one cell to another using this setting. We have tried to reduce the space between the cells by half, but the result was the same. The cell close to the current source would get triggered, but there was no activation of the second cell. See Figure 5.8. Due to the time limitations, we had to stop here with our tests. Although, our current numerical experiments are very limited, GJ coupling seems necessary for the AP propagation. There are some possible explanations that we still would need to test, i.e. smaller mesh step or the distance between the cells. We should note that we have used large distance between cells. We have done this to reduce the computational and time cost of the simulation. But, considering that the GJCs are actually of the size of few nano meters, this means that the distance between the adjacent cells is rather of the size of the cell membrane.

Discussion

Gap junctions play an important role for signal propagation in cardiac tissue. They have a non-linear dynamics that is neglected in the current mathematical models. There are experimental evidences that alternations in the expressions of their connexins affect the propagation velocity. In this project we use the experimental data to fit the non-linear 0D model for homotypic Cx43 and Cx45 GJCs. We propose a discrete spatial model for studying effects of non-linear GJCs in the tissue, but we do not obtain any conclusive results.

Gap junctions are channels of a very small size, and cells almost touch each others in the regions of gap junctions. This implies that our model with the physical spacing might not be appropriate. One idea to modify a model would be to let cells share the boundary, and impose the dynamics of GJCs as a boundary conditions between cells in a same way the ionic models are imposed between the intra-and extracellular space. In any case, there is still a lot of thinking and work to be done on this topic before we can make any conclusion.

. on tissue modelling

One of the open questions in the cardiac electrophysiology consists in understanding the effects of the microscopic heterogeneities of the tissue on the action potential propagation. Standard way to model these is to use the bidomain model and tune the intracellular and extracellular conductivities, based on the intuition.

In this thesis we have proposed a model that extends the standard bidomain model with the periodic diffusive inclusions. These diffusive inclusions represent non-excitable regions in the cardiac tissue, such as fibrotic or fatty tissue. The scale of such diffusive inclusions is significantly larger than the scale of the extracellular space between the neighbouring cells, but much smaller compared to the scale of the tissue.

From a computational point of view this microscopic model is very costly, as we need to crate very detailed meshes of the tissue with the diffusive sub-domains, and the mesh step would depend directly on the size of the diffusive subdomains and the periodic cell. To avoid this problem, we use the technique of multiscale mathematical analysis, the homogenisation theory, to derive an averaged macroscopic model. The homogenised model is in fact the bidomain model, where the effects of the microstructure are observed within the modified intracellular and extracellular conductivity tensors. In the derivation of the modified conductivity tensors we see that both volume and shape of the diffusive inclusions play a role.

To quantify these effects, we have performed 2D and 3D numerical tests. We observed that the modified conductivity tensors, for different shapes, will express different dependence on the volume fractions. Additionally, we observed changes in the anisotropy ratios for both intracellular and extracellular conductivities. Tests on the wave propagation have shown that the inclusions modify the velocity as well as the shape of the wavefront. An interesting result is that for some shapes and volume fractions, the principal direction of the propagation changes. When the inclusions have a shape that is symmetric w.r.t. the fiber direction, the principal directions either stay the same, or change to the perpendicular one. On the other hand, when the inclusions are not aligned with the fiber directions, the principal direction might be under different angle to the fiber direction. Hence, the knowledge of the fiber direction might not be sufficient parameter to determine the main direction of the AP propagation. We conclude that also quantitatively both, volume and geometry of inclusions can play a role in the AP propagation.

To demonstrate that our model is applicable in real-life settings, we simulated the AP propagation in a rat heart, for which enough imaging data were available. As a starting point we used data from HR-MRI of the rat heart, together with data from DT MRI on the fiber directions of the heart. We have segmented a part of the left ventricle of the heart, as the computational domain. We analysed the images to obtain an insight into the fibrotic tissue patches in the heart, and to compute the volume fraction as a smoothly changing parameter throughout the tissue. We have created the mesh of the domain. We used the informations on fiber directions and volume fraction of the fibrotic tissue to compute and map the effective conductivities on the mesh. With this information, we run the simulation of the bidomain model, without additional computational cost, compared to the standard bidomain model.

Finally, we conclude that we have a rigorous and practical way to link the effects of the structural changes in the tissue, as observed in HR MRI, to the changes in the macroscopic electrical behaviour.

... on gap junctional modelling

Gap junctions play an important role for signal propagation in cardiac tissue. They have a non-linear dynamical behaviour that is neglected in the current mathematical models. There are experimental evidences that alternations in the expressions of their connexins affect the propagation velocity. In this project we use experimental data to fit a non-linear 0D model for homotypic Cx43 and Cx45 GJCs, and propose a discrete spatial model for studying effects of non-linear GJCs in the cells strand. The numerical computations of such model are very challenging and could not be completed in a given time.

Outlook

From the modelling point of view there are several limitations of the proposed model. The first one is that the inclusions that we address are purely diffusive. This is a good approximation for tissue patches that are fully made of collagen or fat. For some of the non-excitable regions one can still expect to have different kind of cells inside. In that case, one needs to take into account the cells membranes and different kinds of ionic activity.

Another modelling assumption we used are standard boundary conditions on the boundary between the diffusive inclusions and the bidomain model. These are used in coupling torso and the heart, when the torso is considered to be a diffusive domain. It is not clear that these boundary conditions are appropriate for our case, and different kinds might be discussed.

Additionally, for the sake of mathematical 'simplicity' we have assumed the periodic structure of the tissue. This is not very realistic, and one might explore the non-periodic structures and non-periodic homogenisation techniques.

Of course, also our simulations have limitations. We have tested limited number of shapes for the inclusions, and just a few examples of tissue structure and scar shapes. Already, these were enough to observe complex changes in the wavefront. Obviously, adding more complexity in the shapes of inclusions, or in the regions of the tissue where we apply the modified instead of the standard conductivities, would give more complex patterns.

Some of the first things that can be done is to perform more testing where one might vary: shapes of the inclusions, conductivity in the inclusions, or shapes of the scar-like regions.

What we have done in our simulations on the rat heart is to assume the shape of the inclusions to be everywhere the same. This is a significant simplification in the process. One might perform a detailed image analysis, using HR MRI with a high resolution, to asses the actual or approximate shapes of the inclusions. Such approach would be more accurate, but it would certainly add a lot to the complexity to the pre-processing work. However, once this work would be done, the simulation of the homogenised model would be of the same cost as the one for the standard bidomain model.

Finally, the part that is missing in our work is the validation of the model. To perform this we would need to have data not only from HR MRI and DT MRI of the given heart, but also data on the electrical activity in a controlled environment of the same heart.

More advanced modelling of the fibrotic, fatty or scared tissue could be done with our approach if we had a more precise answers on the following questions:

• What are the typical shapes and sizes of the non-excitable regions?

• Do they usually align with fibers and sheets, or not? Can we have an estimate on the distribution of these parameters? • What is an estimate on the conductivity tensor in these regions? • How does the tissue structure changes from the center of the scar to the border end?

Can we obtain more quantitative details from histology?

3 . 3

 33 Case I, σ d = 0.2[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.76. Plots on the right give the change in anisotropy ratios. . . . 3.4 Case I, σ d = 1.5[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.76. Plots on the right give the change in anisotropy ratios. . . . 3.5 Case I, σ d = 3.0[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.76. Plots on the right give the change in anisotropy ratios. . . . 3.6 Case I, σ d = 6[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.76. Plots on the right give the change in anisotropy ratios. . . . 3.7 Case II, σ d = 0.2[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.39. Plots on the right give the change in anisotropy ratios. . . . 3.8 Case II, σ d = 3.0[0.1S/m]. ThPlots on the left give the values for effective intra-and extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.39. Plots on the right give the change in anisotropy ratios. 3.9 Case III, σ d = 0.2[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.39. Plots on the right give the change in anisotropy ratios. . . . 3.10 Case III, σ d = 3.0[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.39. Plots on the right give the change in anisotropy ratios. . . . 3.11 Case IV, σ d = 0.2[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.98. Plots on the right give the change in anisotropy ratios. . . . 3.12 Case IV, σ d = 3.0[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.98. Plots on the right give the change in anisotropy ratios. . . . 3.13 Case V, σ d = 0.2[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.49. Plots on the right give the change in anisotropy ratios. . . . 3.14 Case V, σ d = 3.0[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.49. Plots on the right give the change in anisotropy ratios. . . . 3.15 Case VI, σ d = 0.2[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.49. Plots on the right give the change in anisotropy ratios. . . . 3.16 Case VI, σ d = 3.0[0.1S/m]. Plots on the left give the values for effective intraand extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.49. Plots on the right give the change in anisotropy ratios. . . . 3.17 Case VII, σ d = 0.2[0.1S/m]. Plots on the left give the eigenvalues for effective intra-and extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.63. Plots on the right give the change in anisotropy ratios. 3.18 Case VII, σ d = 3.0[0.1S/m]. Plots on the left give the eigenvalues for effective intra-and extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.63. Plots on the right give the change in anisotropy ratios. 3.19 Case VII results: the eigenvectors of the conductivity tensors, for σ d = 0.2[0.1S/m]. Black arrows correspond to the eigenvectors for intracellular conductivity tensors, and blue for the extracellular. . . . . . . . . . . . . . . 3.20 Case VII results: the eigenvectors of the conductivity tensors, for σ d = 3.0[0.1S/m]. Black arrows correspond to the eigenvectors for intracellular conductivity tensors, and blue for the extracellular. . . . . . . . . . . . . . . 3.21 Unit cell geometries used to test the effects of the diffusive inclusions on the conductivity tensors in 3D, corresponding to the cases in the Table 3.3. . . . 3.22 3D-results for the effective conductivities and anisotropy ratios. Case I: the diffusive inclusions are superellipsoids with n = 20, r x = r y = r z . And the conductivity σ d kk = 0.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.23 3D-results for the effective conductivities and anisotropy ratios. Case I: the diffusive inclusions are superellipsoids with n = 20, r x = r y = r z . And the conductivity σ d kk = 1.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 3D-results for the effective conductivities and anisotropy ratios. Case I: the diffusive inclusions are superellipsoids with n = 20, r x = r y = r z . And the conductivity σ d kk = 3.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.25 3D-results for the effective conductivities and anisotropy ratios. The diffusive inclusions are superellipsoids with n = 20, and Case II: r x = 0.498, r y = r z , the conductivity σ d kk = 3.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.26 3D-results for the effective conductivities and anisotropy ratios. The diffusive inclusions are superellipsoids with n = 20, and Case III: r x = r y = 0.498, the conductivity σ d kk = 3.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.27 3D-results for the effective conductivities and anisotropy ratios. The diffusive inclusions are superellipsoids with n = 20, and Case IV: r y = r z = 0.498, the conductivity σ d kk = 3.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.28 Circular waves. Initial condition on the transmembrane voltage, v. Gating variable set on h 0 = 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.29 Central stimulus. Standard bidomain model. Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms (left to right). . . . . . . . . . . 3.30 Central stimulus. Case I. σ d l,n = 3.0. |Y B | = 0.23(up) and 0, 61(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms. . . . 3.31 Central stimulus. Case I. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms. . . . . . . . . . . . . . . . . . . 3.32 Central stimulus. Case IV. σ d l,n = 0.2. |Y B | = 0.02(up) and 0.61(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms. . . . 3.33 Central stimulus. Case II. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms, 100ms, 150ms and 200ms. . . . . . . . . . . 3.34 Central stimulus. Case V. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms, 100ms, 150ms and 200ms. . . . . . . . . . . 3.35 Central stimulus. Case III. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms. . . . . . . . . . . . . . . . . . . 3.36 Central stimulus. Case VI. σ d l,n = 0.2. |Y B | = 0.51(up) and |Y B | = 0.61(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms. . . . 3.37 Central stimulus. Case VII. σ d l,n = 0.2. |Y B | = 0.37(up) and 0, 6(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms. . . .

3. 41

 41 Conductivities over the x-axis for the given scars. Top left: conductivities for the SCAR L and the inclusions of type IV. Top right: conductivities for the SCAR M and the inclusions of type IV. Bottom left: conductivities for the SCAR S and the inclusions of type IV. Bottom right: conductivities for the SCAR M and the inclusions of type VI. Legend: purple -σ e l * , orange -σ e n * , blue -σ i l * , green -σ i n * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig. 1 . 1 :

 11 Fig. 1.1: Heart blood circulation system. Created by Eric Pierce, CC BY-SA 3.0, commons.wikimedia.org/w/index.php?curid=830253.

Fig. 1 . 2 :

 12 Fig. 1.2: Electrical conduction system of the heart: 1-Sinoatrial node, 2-Atrioventricular node, 3-Bundle of His, 4-Left bundle branch, 5-Left posterior fascicle, 6-Leftanterior fascicle, 7-Left ventricle, 8-Ventricular septum, 9-Right ventricle, 10-Right bundle branch. By J. Heuser and Patrick J. Lynch, CC BY 2.5,commons.wikimedia.org/w/index.php?curid=1734607.

Fig. 1 . 3 :

 13 Fig. 1.3: Flowchart on the role of Mathematical Modelling in the scientific studies.

Fig. 1 . 4 :

 14 Fig. 1.4: Deaths by cause in women (top) and in men (bottom) in Europe. The latest available year, 2012. From [83].

Fig. 1 . 6 :

 16 Fig. 1.6: Action potential of a ventricular muscle cell. From [29].

Fig. 1 . 7 :

 17 Fig. 1.7: The equivalent electrical circuit for a membrane with three types of ionic channels and a leakage current. Created by Nrets, CC BY-SA 3.0, https://commons.wikimedia.org/w/ index.php?curid=17095070.

Fig. 1 . 8 :

 18 Fig. 1.8: Microscopic structure used to derive the macroscopic, bidomain model.

Fig. 1 . 9 :

 19 Fig. 1.9: Reconstruction of the fiber orientation in the rabbit heart using DT-MRI. (A) Lon-axis cross section showing fiber helix angle. (B) Transmural fiber helix angles extracted from a sector of the left ventricular wall. Taken from [23].

Fig. 1 . 10 :

 110 Fig. 1.10: Details of the gap junction: each channel is made of two hemichannels or connexons.Each connexon is running transverse to the opposing membrane. Connexons of the adjacent cells align themselves to form the structural continuity. Edited image taken from[START_REF] Makowski | Gap junction structures: Analysis of the x-ray diffraction data[END_REF].

Fig. 1 . 11 :

 111 Fig. 1.11: Configurations of connexons and gap junction channels for two different connexins (Cx).Fourteen different connexons can be formed when two connexins are co-expressed in a cell: 2 homomeric connexons, consisting of one type of connexin, and 12 heteromeric connexons consisting of two types of connexins. Gap junction channels formed after docking of two connexons can be homotypic when the connexons are identical (composed of one or two types of connexins) or heterotypic when the connexons are different (composed of one or two types of connexins). Taken from[START_REF] Desplantez | Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties[END_REF].

Fig. 1 . 12 :

 112 Fig. 1.12: Left: Confoncal microscopic projection image of an optical section series through a single isolated ventriculat myocyte labeled for Cx43, showing the positions of clusters of transversely oriented GJCs (arrows). Scale bar = 20 µm. Taken from [88].Right: Lattice of regular hexagons, representing the transverse-section of GJ used in the model in[START_REF] Wilders | Limitations of the dual voltage clamp method in assaying conductance and kinetics of gap junction channels[END_REF]. The diameter of a gap junction is d, the diameter of a gap junction channel is d c ≈ 2 nm, and the lattice constant is d hex ≈ 9 nm.

Fig. 1 . 13 :

 113 Fig. 1.13: Left: DIC image (Differential Interference Contrast) of a pair of cells imaged during dual whole cell patch clamp experiment. Taken from [30]. Right: The equivalent electrical circuit of the dual whole cell patch clamp experiment.

Fig. 2 . 1 :

 21 Fig. 2.1: On the left: the idealised full 2D domain, Ω. On the right: the periodic cell, Y .

Figure 2 . 1 ∞Fig. 2 . 2 :

 2122 Fig. 2.2: Plots for v ε , for different size of ε, and for the homogenised v. We can observe the convergence of the solutions.

  and √ σ d their respective symmetric definite positive square roots.

. 62 )

 62 Now, using (2.58), the exponential part is bounded with,

. 70 )

 70 Definition 10. (m-dissipative operator.) The operator (A, D(A)) is m-dissipative in the Banach space X if (i) u -λAu ≥ u , for all λ > 0 and u ∈ D(A), (ii) for all λ > 0, and all f ∈ X, there exists u ∈ D(A) such that u -λAu = f . Theorem 11. (Hille-Yosida-Phillips Theorem.) A linear operator A is the generator of a contraction semigroup in X if and only if A is m-dissipative with dense domain.

Proposition 24 .

 24 for any v ∈ R. Both inequalities fall directly from this. If the initial function v 0 ε (x) satisfies uniform bounds v 0 ε L 2 (Ω B ε ) ≤ C, and ∇v 0 ε L 2 (Ω B ε ) ≤ C, and if the initial function h 0 ε (x) takes values in (0, 1) for all x ∈ Ω B ε , then we have the uniform bounds:

3 84 3. 3 . 1

 8431 Quantitative effects of diffusive inclusions" "One accurate measurement is worth a thousand expert opinions."-Grace Hopper (Programmer, Mathematician and Admiral) Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.2 Effective conductivities . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.2.1 Settings for 2D cell problems . . . . . . . . . . . . . . . . . . . 65 3.2.2 Results in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.2.3 Settings for 3D cell problems . . . . . . . . . . . . . . . . . . . 79 3.2.4 Results in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.3 Wave propagation in 2D tissue patch . . . . . . . . . . . . . . . . . . . . Settings for homogenised problem . . . . . . . . . . . . . . . . 84 3.3.2 Circular wave propagation . . . . . . . . . . . . . . . . . . . . . 85 3.3.3 Spiral wave dynamics . . . . . . . . . . . . . . . . . . . . . . . . 88 3.3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

Fig. 3 . 1 :

 31 Fig. 3.1: Unit cell geometries, symmetric w.r.t. x-and y-axis, used to test the effects of the diffusive inclusions on the conductivity tensors in 2D. Case I: the shape of the inclusion is a circle. The volume of the inclusion can go from 0% to 76% of the unit cell. Case II and III: the shape of the inclusion is an ellipse with r x : r y = 2 : 1 and 1 : 2, resp. The volume of the inclusion can take up to 39% of the unit cell. Case IV: the shape of the inclusion is a superellipse with e = 0.2 and r x : r y = 1 : 1 . The volume of the inclusion can take up to 98% of the unit cell. Case V and VI: the shape of the inclusion is a superellipse with e = 0.2, and r x : r y = 2 : 1 and 1 : 2, resp . The volume of the inclusion can take up to 49% of the unit cell.

Fig. 3 . 2 :

 32 Fig. 3.2: Unit cell with the nonsymmetric inclusion. The volume of the inclusion can take up to 63% of the unit cell.

  -3.18. We show changes of σ i * and σ e * for each case, depending on the volume fraction of Y B and for different σ d . The tested values for σ d are σ d 11 = σ d 22 = 0.02, 0.15, 0.3 and 0.6 S m -1 , for the Case I. For the other cases we have performed the computations only for 0.02 S m -1 and 0.3 S m -1 . Note that in the computations we scale all conductivities with ς =0.1 S m -1 .

  the Case I the extracellular AR changes from 1.98 for |Y D | = 0, to ≈ 1.0 when |Y D | = 0.76. The intracellular AR goes from 9.0 for |Y D | = 0, reaches minimum ≈ 7.0 at |Y D | = 0.5, and then again grows for the larger values of volume fraction |Y D |. See Figures 3.3

  -

  -3.18 and Figures 3.19-3.20.

Fig. 3 . 3 :Fig. 3 . 4 :Fig. 3 . 5 :Fig. 3 . 6 :Fig. 3 . 7 :Fig. 3 . 8 :Fig. 3 . 9 :Fig. 3 . 10 :Fig. 3 . 11 :Fig. 3 . 12 :Fig. 3 . 13 :Fig. 3 . 14 :Fig. 3 . 15 :Fig. 3 . 16 :

 33343536373839310311312313314315316 Fig. 3.3: Case I, σ d = 0.2[0.1S/m]. Plots on the left give the values for effective intra-and extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.76. Plots on the right give the change in anisotropy ratios.

Fig. 3 . 17 :Fig. 3 . 18 :

 317318 Fig. 3.17: Case VII, σ d = 0.2[0.1S/m]. Plots on the left give the eigenvalues for effective intra-and extracellular conductivities with respect to volume ratio of Y D , ranging from 0.0 -0.63. Plots on the right give the change in anisotropy ratios.

Fig. 3 . 19 :

 319 Fig. 3.19: Case VII results: the eigenvectors of the conductivity tensors, for σ d = 0.2[0.1S/m].Black arrows correspond to the eigenvectors for intracellular conductivity tensors, and blue for the extracellular.

Fig. 3 . 20 :

 320 Fig. 3.20: Case VII results: the eigenvectors of the conductivity tensors, for σ d = 3.0[0.1S/m].Black arrows correspond to the eigenvectors for intracellular conductivity tensors, and blue for the extracellular.

Fig. 3 . 22 :Fig. 3 . 23 : 1 Fig. 3 . 24 :Fig. 3 . 25 :Fig. 3 . 26 :Fig. 3 . 27 :

 3223231324325326327 Fig. 3.22: 3D-results for the effective conductivities and anisotropy ratios. Case I: the diffusive inclusions are superellipsoids with n = 20, r x = r y = r z . And the conductivity σ d kk = 0.2.

Fig. 3 . 28 :

 328 Fig. 3.28: Circular waves. Initial condition on the transmembrane voltage, v. Gating variable set on h 0 = 0.9.

3 . 1 .

 31 All seven cases we test with the volume fraction |Y B | = 0.61 and σ d l,n = 0.2. Additionally, for the Case I we run the simulations for σ d l,n = 3.0 and |Y B | = 0.23 and 0.61, and in the cases IV, VI and VII, we run simulations for σ d l,n = 0.2 and |Y B | = 0.02, 0.51, and 0.37, respectively (these are the minimum values for |Y B | that are possible to obtain with respective shapes of the inclusions). One can notice that in the cases I and IV, (Figures 3.30 -3.32), the volume fraction affects mostly the speed of propagation while the shape of the wavefront stays similar to the one in the standard case. Interestingly, in these two cases σ d does not seem to play an important role (Figure 3.36). It seems that the wavefront for the case I and the volume fraction |Y B | = 0.23, is almost the same as for the case IV and the much smaller volume fraction |Y B | = 0.02.

Figure 3 .

 3 35 and 3.36), one can observe the switch in the principal direction of the propagation. While in the case III we observe it for |Y B | = 0.61, for the case VI we need to reduce |Y B | to 0.51.

Fig. 3 . 29 :

 329 Fig. 3.29: Central stimulus. Standard bidomain model. Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms (left to right).

Fig. 3 . 30 :

 330 Fig. 3.30: Central stimulus. Case I. σ d l,n = 3.0. |Y B | = 0.23(up) and 0, 61(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms.

Fig. 3 . 31 :

 331 Fig. 3.31: Central stimulus. Case I. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms.

Fig. 3 . 32 :

 332 Fig. 3.32: Central stimulus. Case IV. σ d l,n = 0.2. |Y B | = 0.02(up) and 0.61(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms.

Fig. 3 . 33 :

 333 Fig. 3.33: Central stimulus. Case II. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms, 100ms, 150ms and 200ms.

Fig. 3 . 34 :

 334 Fig. 3.34: Central stimulus. Case V. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms, 100ms, 150ms and 200ms.

Fig. 3 . 35 :

 335 Fig. 3.35: Central stimulus. Case III. σ d l,n = 0.2. |Y B | = 0.61. Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms.

Fig. 3 . 36 :

 336 Fig. 3.36: Central stimulus. Case VI. σ d l,n = 0.2. |Y B | = 0.51(up) and |Y B | = 0.61(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms.

Fig. 3 . 37 :

 337 Fig. 3.37: Central stimulus. Case VII. σ d l,n = 0.2. |Y B | = 0.37(up) and 0, 6(down). Observing transmembrane voltage, v at times: 20ms, 50ms and 100ms.

Fig. 3 . 38 :

 338 Fig. 3.38: Initial condition for the simulation where we try to obtain the spiral. The domain size is [-30, 30] × [-[START_REF] Del Corsso | Transfection of mammalian cells with connexins and measurement of voltage sensitivity of their gap junctions[END_REF][START_REF] Del Corsso | Transfection of mammalian cells with connexins and measurement of voltage sensitivity of their gap junctions[END_REF]. On the left is the scaled transmembrane voltage v 0 , which is zero everywhere except in the small region -30 < x < -17 and -3 < y < 3. On the right is the gating variable h 0 that takes a constant value 0 under x-axis, and the value h 0 = 0.5 above x-axis.

Fig. 3 . 39 :

 339 Fig. 3.39: Simulation of the standard bidomain model. The tissue is assumed to be homogeneous, with constant conductivities. We observe the beginning of the spiral, but it quickly disappears. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms, ... , 480ms.

Fig. 3 .

 3 Fig. 3.40: SCARS: The heterogeneous structure of the tissue. We plot the volume fraction |Y B | that ranges from 0.0 for the "dead" center of the scar, to 1.0 for the healthy tissue. The rings where the value ranges 0 < |Y B | < 1, we call the border zone. On the left (SCAR L) the center of the scar has radius r 1 = 16, in the middle (SCAR M) r 1 = 9, on the right (SCAR S) r 1 = 4. For all tree cases the whole scar has the radius r 2 = 22.

Fig. 3 . 41 :.Fig. 3 . 42 :

 341342 Fig. 3.41: Conductivities over the x-axis for the given scars. Top left: conductivities for the SCAR L and the inclusions of type IV. Top right: conductivities for the SCAR M and the inclusions of type IV. Bottom left: conductivities for the SCAR S and the inclusions of type IV. Bottom right: conductivities for the SCAR M and the inclusions of type VI. Legend: purple -σ e l * , orange -σ e n * , blue -σ i l * , green -σ i n

Fig. 3 . 43 :

 343 Fig. 3.43: SCAR M. Case IV. Simulation of the modified bidomain model, with the inclusions assumed to be of the type IV. The tissue is assumed to be heterogeneous as in SCAR M. Like before, the scar is big enough to stop the propagation of AP in the center of the tissue, and forces the wavefront to go around the scar. Hence, we observe again the spiral that is guided by the center of the scar. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 60ms, 120ms, 180ms, ... , 1680ms.

Fig. 3 . 44 :

 344 Fig. 3.44: SCAR S. Case IV. Simulation of the modified bidomain model, with the inclusions assumed to be of the type IV. The tissue is assumed to be heterogeneous as in SCAR S.In this case the center of the scar is small, and not enough to stop or slow down the propagation of AP trough the center. Although, we do observe a change in the shape of the wavefront, we lose the spiral. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms, ... , 480ms.

Fig. 3 . 45 :

 345 Fig. 3.45: SCAR M. Case VI. Simulation of the modified bidomain model, with the inclusions assumed to be of the type VI. The tissue is assumed to be heterogeneous as in SCAR M.There is a change of the direction of the propagation of AP in the scar, which alternates the shape of the wavefront significantly. This kind of activation does not provide the spiral in the tissue in contrast to the same setting with the inclusions of the type IV. Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms, ... , 480ms.

(4. 4 ) 4 . 2

 442 Modified bidomain model adapted for the rat heart

τ

  in = 0.073, τ out = 8.369, τ open = 25.743, τ close = 15.438, v gate = 0.02.

Fig. 4 . 1 :

 41 Fig. 4.1: Flowchart of image analysis.

Fig. 4 . 2 :

 42 Fig. 4.2: HD-MR Images of the rat heart. From Giblert et al. (2012), provided by IHU-Liryc, Bordeaux. The size of the image 256 × 256 × 512. Screenshot of the data imported in Seg3D. Left -volume view; Right -sagital, coronal and axial view of the selected cross-sections.

( a )

 a Segmentation of the slab of the tissue, position in the heart. (b) Segmentation of the slab of the tissue, volume view.

Fig. 4 . 3 :

 43 Fig. 4.3: Segmentation of the tissue slab using Seg3D.

  Figure 4.4.

Fig. 4 . 4 :

 44 Fig. 4.4: Defining inclusions. Pink -the segmented computational domain, purple -inclusions.

Fig. 4 . 5 :

 45 Fig. 4.5: Defining inclusions using different threshold values on the segmentation. Pink -the segmented computational domain. Green -the volume representation of the inclusions.

Fig. 4 . 6 :

 46 Fig. 4.6: Inclusions with the MR Images. Green -the volume representation of the inclusions.

4 .

 4 Assign to each voxel of the image this local volume ratio, |Y D | := num 125 .

Fig. 4 . 7 :

 47 Fig. 4.7: Mesh of the computational domain. Generated with SCIRun. Size: 351706 nodes and 1924747 elements.

Fig. 4 . 8 :

 48 Fig. 4.8: Volume ratio of the inclusions, |Y D |, on the mesh nodes.

Fig. 4 . 9 :

 49 Fig. 4.9: The bidirectional vector field of the fiber orientations. Data obtained from Giblert et al.2012[START_REF] Gilbert | Visualization and quantification of whole rat heart laminar structure using high-spatial resolution contrast-enhanced MRI[END_REF], the original data on the full heart (left) and the projection of the fibers on the segmentation (right). Visualised with SCIRun.

  2.4, we calculated the effective conductivities σ i * l,n and σ e * l,n with respect to the volume fraction, |Y B |, and for given σd . Here we will use σd = 0.2 and 3.0. We already have assigned the local volume fraction, |Y D | = 1 -|Y B |, to each node in the mesh. The plot of |Y B | is shown in the Figure 4.11. Now, we are going to use the results from the section 3.2.4 as a "look-up table" to assign the modified conductivities to each node in the mesh, based on this local volume fraction. Obviously, in the "look-up table" we have the exact results for certain values of |Y B | < 1. For those values of |Y B | for which we do not have the exact results on σ i * l,n and σ e * l,n , we use the linear approximation between two closest known values. Also, for the nodes in which |Y B | = 1, we do not change the conductivity values, i.e. σ i * l,n = σi l,n and σ e * l,n = σe l,n . Resulting modified conductivities are assigned on each node of the mesh, and the results are presented in the the Figures 4.12

  and 4.13.

Fig. 4 . 11 :

 411 Fig. 4.11: Volume fraction |Y B | on the mesh from two sides of the tissue slab.

Fig. 4 . 12 :

 412 Fig. 4.12: Conductivities eigenvalues after applying the volume fraction |Y B | for σ d = 0.2, without effect of fibers.

Fig. 4 . 13 :

 413 Fig. 4.13: Conductivities eigenvalues after applying the volume fraction |Y B | for σ d = 3.0, without effect of fibers.

Figure 4 . 14 .

 414 Then we run two simulations with the modified conductivities, the first one corresponding to σd = 3.0, and the other corresponding to σd = 0.2. The results are shown in the Figures 4.

15

 15

  and 4.16, respectively.

Fig. 4 . 14 :

 414 Fig. 4.14: Simulation of the bidomain model on the slab of tissue. Times: 5, 10, 15, 20, 25, 30, 35, 40, 60 and 80 ms.

Fig. 4 . 15 :

 415 Fig. 4.15: Simulation of the modified bidomain model on the slab of tissue, corresponding to σd = 3.0. Times: 5, 10, 15, 20, 25, 30, 35, 40, 60 and 80 ms.

Fig. 4 . 16 :

 416 Fig. 4.16: Simulation of the modified bidomain model on the slab of tissue, corresponding to σd = 0.2. Times: 5, 10, 15, 20, 25, 30, 35, 40, 60 and 80 ms.

  and 5.2). The data are obtained from the work of Desplantez et al 2004, 2007 and 2011.

Fig. 5 . 1 :

 51 Fig. 5.1: Dependence of the normalised gap junction conductance g j,∞ (lower curves) and g j,0 (upper curves) on junctional potential Vj for homotypic Cx43 and Cx45 channels.

Fig. 5 . 2 :Tab. 5 . 1 :

 5251 Fig. 5.2: Kinetics of inactivation, τ j,∞ as a function of V j for homotypic Cx43 and Cx45 channels.

Fig. 5 Fig. 5 . 4 :Fig. 5 . 5 : 5 . 5 0D simulation 123 5 . 6

 554555556 Fig. 5.3: Experimental results

Ωe

  ∆tσ e ∇u e ∇v e + Γ ion u e v e -Γ ion(u i -V m + ∆tI ion (V m ), h)v e = 0, Ωe ∆tσ i ∇u i ∇v i + Γ ion u i v i -Γ ion

15 :: u 1 e = u e and u 1 i = u i 17 :

 151117 7: Define: A(u e , v e ) := Ωe ∆tσ e ∇u e ∇v e + Γ ion u e v e -Γ ion(u i -V m + ∆tI ion (V m ), h)v e . A(u i , v i ) := Ωe ∆tσ i ∇u i ∇v i + Γ ion u i v i -Γ ion (u e + V m -∆tI ion (V m ), h)v i .8: for t = 0 to T , t+ = ∆t do 9: Solve A(u e , v e ) = 0 10: Solve A(u i , v i ) = 0 11: u 1 e = u e and u 1 i = u i 12: Solve A(u e , v e ) = 0 13: Solve A(u i , v i ) = 0 14: Compute: D i = max |u i -u 1 i | and D e = max |u e -u 1 e | while D i > threshold or D e > threshold do 16Solve A(u e , v e ) = 0 18: Solve A(u i , v i ) = 0 19: Compute: D i = max |u i -u 1 i | and D e = max |u e -u 1 e | 20:

(Fig. 5 . 7 :Fig. 5 . 8 :

 5758 Fig. 5.7: Preliminary numerical results on 4 myocytes.

6 Conclusions"

 6 "One never notices what has been done, one can only see what remains to be done." -Marie Currie (Physicist and Chemist) Contents 6.1 Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.1.1 ... on tissue modelling . . . . . . . . . . . . . . . . . . . . . . . 131 6.1.2 ... on gap junctional modelling . . . . . . . . . . . . . . . . . . 132 6.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.1 Thesis overview 6.1.1 ..

  Bundle of His, 4-Left bundle branch, 5-Left posterior fascicle, 6-Leftanterior fascicle, 7-Left ventricle, 8-Ventricular septum, 9-Right ventricle, 10-Right bundle branch. By J. Heuser and Patrick J. Lynch, CC BY 2.5,commons.wikimedia.org/w/index.php?curid=1734607. . . . . . . . . . . . . . Flowchart on the role of Mathematical Modelling in the scientific studies. . .

	1.1	Heart blood circulation system. Created by Eric Pierce, CC BY-SA 3.0, com-
		mons.wikimedia.org/w/index.php?curid=830253. . . . . . . . . . . . . . .
	1.2	Electrical conduction system of the heart: 1-Sinoatrial node, 2-Atrioventricular
	node, 3-1.3

1.4 Deaths by cause in women (top) and in men (bottom) in Europe. The latest available year, 2012. From [83]. . . . . . . . . . . . . . . . . . . . . . . . .

  Plots to be read row by row from left to right. Observing transmembrane voltage, v at times: 40ms, 80ms, 120ms,

... , 480ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Flowchart of image analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 HD-MR Images of the rat heart. From Giblert et al. (2012), provided by IHU-Liryc, Bordeaux. The size of the image 256 × 256 × 512. Screenshot of the data imported in Seg3D. Left -volume view; Right -sagital, coronal and axial view of the selected cross-sections. . . . . . . . . . . . . . . . . . . . . 4.3 Segmentation of the tissue slab using Seg3D. . . . . . . . . . . . . . . . . . . 4.4 Defining inclusions. Pink -the segmented computational domain, purpleinclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Defining inclusions using different threshold values on the segmentation. Pink -the segmented computational domain. Green -the volume representation of the inclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Inclusions with the MR Images. Green -the volume representation of the inclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Mesh of the computational domain. Generated with SCIRun. Size: 351706 nodes and 1924747 elements. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Volume ratio of the inclusions, |Y D |, on the mesh nodes. . . . . . . . . . . .

  .13 Conductivities eigenvalues after applying the volume fraction |Y B | for σ d = 3.0, without effect of fibers. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Simulation of the bidomain model on the slab of tissue. Times: 5, 10, 15, 20, 25, 30, 35, 40, 60 and 80 ms. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Simulation of the modified bidomain model on the slab of tissue, corresponding to σd = 3.0. Times: 5, 10, 15, 20, 25, 30, 35, 40, 60 and 80 ms. . . . . . 4.16 Simulation of the modified bidomain model on the slab of tissue, corresponding to σd = 0.2. Times: 5, 10, 15, 20, 25, 30, 35, 40, 60 and 80 ms. . . . . . Dependence of the normalised gap junction conductance g j,∞ (lower curves) and g j,0 (upper curves) on junctional potential Vj for homotypic Cx43 and Cx45 channels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Kinetics of inactivation, τ j,∞ as a function of V j for homotypic Cx43 and Cx45 channels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	5.1

. . . . . . . . . . . . . . . . . . . . 4.11 Volume fraction |Y B | on the mesh from two sides of the tissue slab. . . . . . 4.12 Conductivities eigenvalues after applying the volume fraction |Y B | for σ d = 0.2, without effect of fibers. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Transjunctional current simulated with the parameters obtained from the experiments and the data fitting. . . . . . . . . . . . . . . . . . . . . . . . .

  Parameters of the bidomain model for a rat heart. . . . . . . . . . . . . . . 98 5.1 Parameters obtained from the experiments and the data fitting. The units in the table are the following: G j,max in nS, V j,0 in mV , A and b in mV -1 and a in ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

	4.3
	Introduction to Modelling in Cardiac
	Electrophysiology
	" "Mathematicians are a sort of Frenchmen. If you talk to them, they translate it into their own
	language, and then it is immediately something
	quite different."
	-Johann Wolfgang von Goethe
	(Writer)
	xxiii

4.1

Morphometric data in rat ventricular myocytes . . . . . . . . . . . . . . . .

4.2

Action potential characteristics in epicardial and endocardial models, from

[START_REF] Sandeep | A mathematical model of action potential heterogeneity in adult rat left ventricular myocytes[END_REF]

. Contents 1.1 Heart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Physiology and function . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Electrophysiology . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Pathology statistics . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Cell scale modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Cell membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Ion flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Gating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.4 Ionic models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.5 The Mitchell-Schaeffer model . . . . . . . . . . . . . . . . . . . 1.3 Tissue scale modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Bidomain model . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Fibers and anisotropy . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Asymptotic ionic model . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Limitations of the bidomain model and structural heterogeneities 1.4 The gap junctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Gap junctional structure . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Electrical function . . . . . . . . . . . . . . . . . . . . . . . . . .

  .44) where τ close and τ open are the time constants with which the gate closes and opens, and v gate is the threshold value of the transmembrane potential that indicates the change in the dynamics of the gate. The parameters for the ventricular cells proposed in the original paper are: τ in = 0.3 ms, τ out = 6 ms, τ open = 120 ms, τ close = 160 ms, and v gate = 0.13.

1 :

 1 Define meshes on Ω and Ω B ε . 2: Define T and N = T ∆t . 3: Define var. form. spaces H 1 (Ω), H 1 (Ω B ε ) with P1 elements.

4: Declare u i εn , v εn , h εn , φ i , φ ∈ H 1 (Ω B ε

1 :

 1 Define meshes on Y , and Y B . 2: Define var. form. spaces H 1 (Y ), H 1 (Y B ) with P1 elements, both periodic on the outer boundary ∂Y .

3: Declare w i k , φ i ∈ H 1 (Y B ) and w e k , φ i ∈ H 1 (Y ), for k = 1, 2, or k = 1,

2, 3. 4: Solve for w i k :

  = 0.1S m -1 . The reason for this choice of scaling is given in the Section 1.3.3. We test the following scaled values for σ d = {0.2, 1.5, 3.0, 6.0}. We have to solve four cell problems. The first two on Y B ,

	741 0	0 0.1934	and σ e =	3.904 0	0 1.970	,	(3.3)
	that are scaled with ς Y B

  .2.

				σ d 11,22 [S m -1 ]
	Case I 1.0	1 : 1	0.76	0.02, 0.15, 0.3, 0.6
	Case II 1.0 1 : 0.5	0.39	0.02, 0.3
	Case III 1.0	1 : 2	0.39	0.02, 0.3
	Case IV 0.2	1 : 1	0.98	0.02, 0.3
	Case V 0.2 1 : 0.5	0.49	0.02, 0.3
	Case VI 0.2	1 : 2	0.49	0.02, 0.3

Tab. 3.2:

The test cases for geometries of inclusions in 2D. e r x : r y max |Y D |

  rest (mV ) dV /dt max (V /s) PO (mV ) APD 50 (ms) APD 90 (ms) rest -resting membrane potential, dV /dt max -maximum upstroke velocity, PO -peak overshoot, APD 50 -action potential duration (50% repolarisation), APD 90 -action potential duration (90% repolarisation).

	Epicardial	-80.44	145.41	35.75	12.99	39.68
	Endocardial	-80.37	181.41	47.08	37.83	76.43
	V

Table 4 .

 4 3. Parameters of the bidomain model for a rat heart.

	Tab. 4.3: Parameter	Value	units
	V max	0.035	V
	V rest	-0.08	V
	V m	0.125	V
	χ	8 × 10 5	m -1
	C m	0.01	F/m 2
	ς	0.1	S/m
	σ i l σ i t,n	1.741 0.1934	--
	σ e l σ e t,n	3.906 1.970	--
	T	1 • 10 -3	s
	L	1 • 10 -3	m
	N	0.0125	-
	4.2.3 Main Contribution: image-based computation of the modified
	conductivities

  table are the following: G j,max in nS, V j,0 in mV , A and b in mV -1 and a in ms.

		G j,max g j,min V j,0	A	a	b
	Cx43	68.3	0.27 60.5 0.098 33880 -0.06
	Cx45	2.0	0.16 38.6	0.1	132330 -0.07

  Thoughts on the tissue model with non-linear GJCs and we keep the homogeneous Neumann boundary conditions on the rest of the external boundary. We use the following variational formulation of the problem:Ω i σ i ∇u i ∇φ i + Ωe σ e ∇u e ∇φ e -Γ ion (∂ t V m -I ion (V m , h))(φ i -φ e ) -

	Γ j

5.6

1.1 Heart

1.2 Cell scale modelling

1.3 Tissue scale modelling

1.4 The gap junctions

2.2 Microscopic problem statement

2.3 Formal derivation of the macroscopic problem

2.4 The rigorous mathematical proofs: existence and convergence

2.5 Discussion

3.2 Effective conductivities

3.3 Wave propagation in 2D tissue patch

Chapter 3 Quantitative effects of diffusive inclusions

https://www.ihu-liryc.fr/en/ 4.3 Image analysis

4.3 Image analysis

http://teem.sourceforge.net/nrrd/index.html4.3 Image analysis

http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf

http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format

4.4 Simulation in FreeFem++

5.4 Data fitting

-σ i ∆u i =0, in Ω i , -σ e ∆u e =0,in Ω e ,

5.[START_REF] Amar | Existence and uniqueness for an elliptic problem with evolution arising in electrodynamics[END_REF] Thoughts on the tissue model with non-linear GJCs

5.7 Discussion

6.2 Outlook

This thesis was mainly funded by Inria PhD programme, and performed in Carmen team at Inria Bordeaux Sud Ouest Research Center.

This thesis was partly performed in the ihu-Liryc institute for rhythmology and cardiac modelling, and received financial support from the French Government as part of the "Investments of the Future" program managed by the National Research Agency (ANR), Grant reference ANR-10-IAHU-04.

This work was partially funded by ANR HR-CEM project (ANR-13-MONU-0004): for my PhD salary (3 months), for study visit to University Ottawa (Canada), for IBBM summer school at University of Utah (USA), for participation in LACM conference at TUM (Germany).

Funding

the walls of the unit cell, and we vary r y = r z over the range (0, 0.49). In this way the inclusions are elongated, and aligned with fibers. In the last two cases the inclusions have a shape of sheets. In the Case III they are perpendicular to the z-axis, and in the Case IV they are perpendicular to the x-axis.

Tab. 3.3:

The test cases for the geometries of the inclusions in 3D.

Results in 3D

The results are shown in the Figures ?? -3.27. Like before we plot σ i * and σ e * depending on the volume fraction |Y D |. In all of the tested cases the volume fraction can range from 0.02 to 1. For the first case we assume the values for σ d kk = 0.2, 1.5 and 3.0, also scaled with ς = 0.1S/m, and for the other tree cases we choose only σ d kk = 3.0. The results in 3D show the similar behaviour as the ones in 2D. The limit cases are the same as before: when the inclusions are small, i.e. when their volume is close to zero, the effective conductivities are almost the same as the conductivities in the standard bidomain model; on the contrary side when the inclusions are big, i.e when their volume is almost one, the effective conductivities are like in 2D case σ i * is approaching the 0-tensor and σ e * is approaching σ d . This is observed for the Case I, where we have assumed that

Mesh

The next step is to obtain the mesh of the segmented domain. The SCIRun software is organised in modules, where each module performs one task. The modules can be connected to each others, so that the output of one module is used as an input of the other. The programs are essentially schemes, visual and intuitive. The available modules go from those for importing files, via those for performing many functions from simple arithmetic to solving partial differential equations, to those for the visualisation. We will make use of this software for two main things: generating the mesh and mapping the relevant parameters on the mesh. For the mesh generation the main modules that we use are:

• Read Field -module that, as the name says, reads the input file. In this case the segmentation file, generated with Seg3D.

• ExtractIsosurface and FairMesh -modules for setting the parameters for mesh generation. In particular, ExtractIsosurface extracts an isopotential surface from a scalar field. And the FairMesh smooths surface meshes without shrinking them. These two are used for defining the surface that will be filled with tetrahedra with the following module.

• InterfaceWithTetgen -module for calling the TetGen software for meshing, where we can set additional options. Among others we can set the mesh quality by setting the minimum radius-edge ratio, and impose the maximum volume constraint on all tetrahedra.

The generated mesh can be exported with different file extensions, we chose .VTK. The result is a fine mesh, that has: 351706 nodes and 1924747 elements, i.e. tetrahedra. See Figure 4.7.

Mapping volume fraction on mesh

For the numerical simulation, every parameter of the simulation has to be defined on the mesh. At the moment we have the volume fraction defined on each voxel of the image, and this parameter needs to be mapped on the mesh. For this purpose we use again SCIRun software. The modules used here are the following:

• Two different Read Field modules, one reads the mesh file and the other the file with the information about the volume fraction.

• MapFieldDataOntoNodes -Interpolation of data from one mesh or point cloud to another mesh or point cloud. The output mesh will have the data located at the nodes. We use it for mapping the volume fractions defined on the voxels of the

Introduction

Here we deal with the propagation of the signal on the scale of a few cardiac cells. Such models are numerically very challenging, and they should not be used for the simulations on the larger scales. But they provide closer insight into the contributions of different parameters of the models. More specifically, we are interested in understanding the non-static electrical properties of gap junctions.

The cardiac cells are surrounded by an extracellular matrix and neighbouring cells are electrically connected via gap junctions (GJ). Gap junctions are clusters of gap junction channels (GJCs) that are mainly localised on the longitudinal ends of the cells, where they compose the intercalated disks. These channels play a role of resistance pathways for the direct cell to cell propagation by providing the direct passage of molecules and ions between adjacent cells. See Section 1.4.

In the usual continuous models for the electrical propagation, the GJCs are tiny channels between cells which interior domains extend the intra-cellular domain, e.g. simple linear electrical conductors. However, the voltage clamp experiments (e.g. [START_REF] Harris | Kinetic properties of a voltage-dependent junctional conductance[END_REF][START_REF] Del Corsso | Transfection of mammalian cells with connexins and measurement of voltage sensitivity of their gap junctions[END_REF][START_REF] Desplantez | Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties[END_REF]) have shown that GJCs exhibit a non-linear voltage and time dependence on the The equivalent ODE to (5.4), now can be written as

(5.7)

From the kinetic theory and Boltzmann relation, we know that the opening and closing rates, α and β are exponentially related to the energy difference between these two states, i.e. to the voltage field across the junction, and can be expressed as

where V j is the transjunctional voltage, and V j,0 is the transjunctional voltage at which the conductance is half minimal (n = 1/2), and A α and A β are constants expressing the voltage sensitivity. From here and from (5.6), we express n ∞ , as

where A = A α + A β . Now, if γ o and γ c are the conductivities of a single open and closed channel, and if N is the total number of channels, then the maximal conductance is G j,max = γ o N , the minimal conductance is G j,min = γ c N , and total conductance of the gap junction is given as

(5.10)

= n(G j,max -G j,min ) + G j,min .

(5.11)

To obtain the expression for the gating variable we need to scale the total conductance with its maximal value, g j = G j /G j,max . So, the gating variable is given as, g j = n(1 -g j,min ) + g j,min .

(5.12)

The steady state value of the gating variable is g j,∞ = n j,∞ (1 -g j,min ) + g j,min , (5.13) and using (5.9), we have, g j,∞ (V j ) = 1 -g j,min 1 + e A(V j -V j,0 ) + g j,min .

(5.14)

From (5.6),(5.8) and , we can see that the time constant also has an exponential dependence on V j , that we simplify into τ j,∞ (V j ) = ae -bV j .

(5.15)

Colophon

This thesis was typeset with L A T E X 2 ε . It uses the Clean Thesis style developed by Ricardo Langner. The design of the Clean Thesis style is inspired by user guide documents from Apple Inc. Download the Clean Thesis style at http://cleanthesis.der-ric.de/.