Hubert Cochet

Michel Haïssaguerre

Pierre Jaïs

Benjamin Berte

Yuki Komatsu

Rémi Dubois

Erin Stretton

Krissy Macleod

Marine Breuilly

Radiofrequency Ablation Planning for Cardiac Arrhythmia Treatment using Modeling and Machine Learning Approaches

Keywords: Cardiac electrophysiology modeling, intracardiac electrogram modeling, machine learning, radiofrequency ablation planning, electroanatomical mapping

Cardiac arrhythmias are heart rhythm disruptions which can lead to sudden cardiac death. They require a deeper understanding for appropriate treatment planning.

In this thesis, we integrate personalized structural and functional data into a 3D tetrahedral mesh of the biventricular myocardium. Next, the Mitchell-Schaeffer (MS) simplified biophysical model is used to study the spatial heterogeneity of electrophysiological (EP) tissue properties and their role in arrhythmogenesis.

Radiofrequency ablation (RFA) with the elimination of local abnormal ventricular activities (LAVA) has recently arisen as a potentially curative treatment for ventricular tachycardia but the EP studies required to locate LAVA are lengthy and invasive.

LAVA are commonly found within the heterogeneous scar, which can be imaged non-invasively with 3D delayed enhanced magnetic resonance imaging (DE-MRI). We evaluate the use of advanced image features in a random forest machine learning framework to identify areas of LAVA-inducing tissue. Furthermore, we detail the dataset's inherent error sources and their formal integration in the training process.

Finally, we construct MRI-based structural patient-specific heart models and couple them with the MS model. We model a recording catheter using a dipole approach and generate distinct normal and LAVA-like electrograms at locations where they have been found in clinics. This enriches our predictions of the locations of LAVA-inducing tissue obtained through image-based learning. Confidence maps can be generated and analyzed prior to RFA to guide the intervention. These contributions have led to promising results and proofs of concepts.

Planification de l'ablation radiofréquence des arythmies cardiaques en combinant modélisation et apprentissage automatique Résumé:

Les arythmies sont des perturbations du rythme cardiaque qui peuvent entrainer la mort subite et requièrent une meilleure compréhension pour planifier leur traitement.

Dans cette thèse, nous intégrons des données structurelles et fonctionnelles à un maillage 3D tétraédrique biventriculaire. Le modèle biophysique simplifié de Mitchell-Schaeffer (MS) est utilisé pour étudier l'hétérogénéité des propriétés électrophysiologiques (EP) du tissu et leur rôle sur l'arythmogénèse.

L'ablation par radiofréquence (ARF) en éliminant les activités ventriculaires anormales locales (LAVA) est un traitement potentiellement curatif pour la tachycardie ventriculaire, mais les études EP requises pour localiser les LAVA sont longues et invasives.

Les LAVA se trouvent autour de cicatrices hétérogènes qui peuvent être imagées de façon non-invasive par IRM à rehaussement tardif. Nous utilisons des caractéristiques d'image dans un contexte d'apprentissage automatique avec des forêts aléatoires pour identifier des aires de tissu qui induisent des LAVA. Nous détaillons les sources d'erreur inhérentes aux données et leur intégration dans le processus d'apprentissage.

Finalement, nous couplons le modèle MS avec des géométries du coeur spécifiques aux patients et nous modélisons le cathéter avec une approche par un dipôle pour générer des électrogrammes normaux et des LAVA aux endroits où ils ont été localisés en clinique. Cela améliore la prédiction de localisation du tissu induisant des LAVA obtenue par apprentissage sur l'image. Des cartes de confiance sont générées et peuvent êtres utilisées avant une ARF pour guider l'intervention. Les contributions de cette thèse ont conduit à des résultats et des preuves de concepts prometteurs.

Mots clés: modélisation de l'électrophysiologie cardiaque, modélisation d'électrogrammes intracardiaques, apprentissage automatique, planification d'ablation par radiofréquence, activités ventriculaires anormales locales (LAVA)

Clinical Context

Cardiovascular diseases (CVD) refer to the conditions that affect the heart and the circulatory system. To date, they remain the leading cause of death in the western world. According to the Global Burden of Disease, CVD were responsible for more than 29% of deaths in the world in 2013 (> 15,616 million deaths), twice the amount of deaths caused by cancer in the same year [START_REF] Nichols | [END_REF]]. In the European continent alone, despite the recent decrease in CVD deaths, these diseases claim over 4 million victims per year. Figure 1.1 shows the proportion of non-accidental deaths in Europe in 2013. Cardiac arrhythmias are a subset of CVD grouping abnormalities in the heart rhythm. A dangerous consequence of these rhythm perturbations includes a compromise of the heart's effectiveness to pump blood. Sudden cardiac death (SCD) occurs if the condition is not treated within a short delay. Early detection and accurate prediction of disease progression of CVD remain an important need to reduce their mortality. Furthermore, improvements in therapy planning and guidance are of vital importance to reduce the mortality of these diseases.

Depending on the nature of the arrhythmia, either implantable cardioverterdefibrillator (ICD) or radiofrequency ablation (RFA) treatment can be recommended. An ICD is a device that continuously monitors the electrical activity in the heart and delivers electrical shocks to control life-threatening arrhythmias. Nonetheless, their effectiveness in reducing mortality is still unclear [START_REF] Nattel | Decade in reviewarrhythmias: Cardiac fibrillation challenges and evolving solutions[END_REF]]. Besides, this therapy is a non-curative option and is commonly paired with pharmaceutical treatment. A potentially curative alternative is RFA treatment, where thermal lesions are generated in the heart to interrupt abnormal re-entry circuits that cause arrhythmias. The greatest challenge in this therapy is ablation target identification. Currently, this can be achieved using electrophysiological (EP) substrate mapping. Recently, [START_REF] Jaïs | [END_REF] showed that the elimination of local abnormal ventricular activities (LAVA) was an effective endpoint for substrate-based ventricular tachycardia.

Figure 1.1: Proportion of all deaths due to major causes in Europe (Top) among men and (Bottom) among women. Courtesy of [START_REF] Nichols | [END_REF] The main questions we aim to answer in this thesis are:

-Can personalized structural and functional data help us understand the patient-specific arrhythmogenesis?

-Can appropriate imaging data provide us with information to characterize and identify RFA targets?

-Are we able to use non-invasive image-based models and in silico simulation to generate intracardiac electrograms with LAVA-like patterns at the locations where they have been found in clinics?

-Can these simulated electrograms enrich the predictions of RFA targets obtained through image-based learning?

1.2. Manuscript Organization 3

Manuscript Organization

This thesis is based on our published and submitted work, Figure 1.2 exemplifies the manuscript organization. In Chapter 2, we introduce the relevant clinical concepts for understanding this work and review the state of the art of electrophysiology modeling and the use of machine learning schemes in the medical domain.

In Chapter 3, based on [Cabrera-Lozoya 2015c], we make use of integrated personalized structural and functional data to study the spatial heterogeneity of electrophysiological tissue properties and their role in arrhythmogenesis.

In Chapter 4, based on [Cabrera-Lozoya 2014] and on [Cabrera-Lozoya 2015d], we evaluate the predictive power of locally computed intensity and texture-based MRI features to identify radiofrequency ablation targets using a machine learning framework, including a detailed analysis of the dataset's inherent error sources and their integration in the training process.

In Chapter 5, based on [Cabrera-Lozoya 2015a], we test the feasibility of using delayed enhanced MR image-based models to reproduce LAVA-like patterns in simulations of catheter recordings of intracardiac electrograms.

In this thesis, we use mathematical modeling and machine learning techniques alongside with patient-specific clinical data to guide radiofrequency ablation therapy for the treatment of cardiac arrhythmias. In this chapter, we review the essentials and basics of the clinical aspects and mathematical techniques.

The Heart: Generalities

Cardiac Anatomy

The heart is the muscular organ responsible for pumping blood through the cardiovascular system, providing the body with oxygen and nutrients while removing metabolic waste. It is enclosed in a protective sac called the pericardium, which prevents overfilling of the heart and helps it anchor to the surrounding structures in the mediastinum in the chest. The heart is composed of four chambers: the left and right atria (upper, bloodreceiving chambers) and the left and right ventricles (lower, blood-discharging chambers) as shown in Figure 2.1a. Blood from systemic circulation flows through the superior and inferior vena cavae into the right atrium, across the tricuspid valve and into the right ventricle. Then, the blood crosses the pulmonary valve and passes into the pulmonary circulation in direction of the lungs. Once metabolic wastes are removed from the blood and oxygen is added, it will return to the left atrium, through the mitral valve and into the left ventricle before being discharged again to the systemic circulation through the aorta.

The heart wall is divided in three layers shown in Figure 2.1b, the outer-most layer is called the epicardium, the middle one composed mainly of cardiac muscle responsible for heart contraction is referred to as the myocardium and the inner layer, the endocardium, is in direct contact with the blood. Heart muscle cells, myocytes, are arranged in fibers whose direction varies across the myocardium. The fiber orientation can be described by an elevation angle corresponding to its obliquity with respect to the plane of the section. The value of this angle varies from approximately +70 • on the endocardium, 0 • in mid-myocardium to -70 • in the epicardium (Figure 2.2). This swirled organization optimizes the pumping function of the heart.

Cardiac Electrophysiology in the Healthy Heart

An adult human heart at rest contracts with a stable rhythm of approximately 60-100 beats per minute (bpm). This sinus rhythm is managed by the automaticity of the heart's natural pacemaker, the sinoatrial (SA) node, which is a group of cells capable of spontaneously generating electrical signals that propagate through the heart muscle. The electrical wave travels first through the atrial fibers (the pathway propagating the electrical impulse to the left atrium is known as the Bachmann's bundle) and to the atrioventricular (AV) node. Here, a delay allowing an appropriate filling of the ventricles during the atrial contraction takes place. Thereafter, the electrical stimulus travels through the His bundle into the left and right bundle branches. Finally, the Purkinje network allows the activation of the cardiac myocytes in the ventricles. A schematic of the heart's conduction system is found in Figure 2.3. In non-pathological conditions, a cardiac cell has a resting polarized state, meaning that there exists an electrical potential difference, typically of -90mV, across the cell membrane separating the intracellular and extracelluar mediums. When the cell receives an electrical stimulus (like the wave arising from the SA node), its transmembrane potential evolves in accordance with the gradient of ionic concentrations inside and outside the cell. The temporal evolution of the transmembrane potential Chapter 2. Clinical and Technical Context is known as an action potential and it consists of three main phases (Figure 2.4):

1. Depolarisation: Fast N a + channels open and sodium ions can enter the cell, rapidly rendering its inside positive with respect to the extracellular medium.

2. Plateau: K + and N a + channels are slowly de-activated and Ca 2+ ones are opened, momentarily stabilizing the transmembrane potential.

3. Repolarisation: The K + channels gradually re-activate and there is a decrease in the calcium influx, re-establishing the ionic rest state of the cell.

The Pathological Heart

Myocardial Infarctions

In pathological cases, an interruption of the blood supply to the cardiac tissue can occur. The main cause of this phenomenon is the occlusion of one of the blood vessels irrigating the heart, known as the coronary arteries. If left untreated, oxygen-deprived cells can suffer from irreversible damage or even death (myocardial infarction, MI). Areas adjacent to the necrotic core are referred to as the border zone, where surviving myocytes are surrounded by dense fibrosis, as seen in Figure 2.5. These areas present decreased lateral connection between the myocardial cells, which undergo changes in the transmembrane ion channel expression leading to differences in their electrophysiological behaviors [START_REF][END_REF]]. These disruptions can potentially cause abnormal heart beating patterns (cardiac arrhythmias) or even a complete stopping of the heartbeat, a condition known as cardiac arrest.

Cardiac Arrhythmias

A cardiac arrhythmia is a broad classification of the conditions that perturb the normal cardiac rhythm, rendering it either irregular, too slow or too fast. Of particular interest for this thesis are ventricular tachycardias, which are characterized by an abnormally fast (from 120 up to 280 bpm) heart rhythm originating

Scar-related ventricular tachycardia (VT) or arrhythmia

The presence of scar in the myocardial tissue can lead to abnormal electrical circuit generations, often referred to as scar-related reentry [Aliot 2009]. Although the main cause for scar formation reported is myocardial infarction [Stevenson 1993], other conditions can result in damaged myocardial tissue, including arrhythmogenic right ventricular cardiomyopathy (ARVC) [Belhassen 1984], dilated cardiomyopathy [Hsia 2002] or surgery [Eckart 2007]. Myocardial scar can be identified due to its electrophysiological characteristics, which include low-voltage regions or fractionated electrograms, or through imaging studies. According to [Aliot 2009], substrate promoting scar-related re-entry can have the following characteristics:

-Regions of slow conduction.

-Unidirectional conduction block at some point in the re-entry path that allows initiation of re-entry.

-Areas of conduction block that often define parts of the re-entry path.

According to animal [El-Sherif 1981] [El-Sherif 1983] and human [De Bakker 1988] studies, the phenomenon of re-entry is due to the existence of surviving myocyte bundles and the increased presence of collagen and connective tissue affecting their coupling with the rest of the myocardium. This perturbation of normal myocardial architecture promotes abnormal patterns of activation and conduction velocity [Peters 1998].

Imaging the Heart

In patients with suspected scar-related VT, imaging studies can be used to assess the viability of the myocardial tissue, as well as to obtain important aspects of tissue heterogeneity, scar extent and location. Delayed-enhanced (DE) magnetic resonance imaging (MRI), which has the advantage of being a non-ionizing imaging technique, is used in conjunction with gadolinium (Gd) as a contrast agent to improve scar visualization. Nonetheless, DE-MRI is often prohibited in VT-prone patients because many of them have implanted pacemakers or ICDs. Also, the presence of these devices can lead to important artifacts and deteriorate image quality. The use of DE multi-detector computed tomography (MDCT) has been investigated in these cases, but its reproducibility is still debated [Aliot 2009].

To date, DE-MRI remains the gold standard for myocardial size and morphology evaluation, it is also used to assess the heterogeneity of the border zone [Schuleri 2009]. Figure 2.6 exemplifies the use of DE-MDCT and DE-MRI in animal subjects for visualizing the scar and border zone (also known as peri-infarct zone) and their correlation with histological cuts.

Electro-anatomical mapping

Electroanatomical mapping (EAM) is a minimally-invasive technique used to record in-vivo cardiac electrical activity at specific locations inside the heart. The method of entry is chosen according to the cardiac chamber to be studied. When the cavity of study is the right part of the heart, then the mapping catheter is introduced through the femoral vein. Alternatively, when studying the left heart, the femoral artery is chosen. When the catheter either comes in contact with the tissue of interest or is located within its vicinity, both the electrical activity of the tissue and the coordinates of the catheter's position in space are retrieved. These systems integrate three main functions [Aliot 2009 Despite the capability of locating the catheter in 3D space and therefore reconstructing the heart chamber that is being mapped, geometries provided by EAM systems tend to be rough estimates of the actual cardiac anatomy. Catheter position is highly affected by cardiac or respiratory motion, and anatomic reconstruction algorithms may vary between systems [Aliot 2009]. These limitations can compromise the integration of EAM information with anatomical information from traditional imaging systems. Often, multiple mapping catheters are used in order to improve spatial sampling.

Mapping Systems

Two main types of mapping systems currently exist in the clinical market: contact mapping and non-contact mapping. Contact Mapping: An example of this type of mapping system is the CARTO As the name implies, this mapping technique requires direct contact between the recording catheter and the cardiac tissue. The catheter is moved along the endocardial or epicardial wall to capture electrical signals. It is preferred for the study of local phenomena. The three-dimensional location of the catheter in this particular system is achieved through the placement of three coils beneath the patient producing low-level electromagnetic fields which are measured by a sensor located at the tip of the mapping catheter, which have reported a magnitude of location errors from the catheters of ∼ 1mm [Gepstein 1997]. In some applications throughout this thesis, the contact mappping system was used in conjunction with a multi-spline recording catheter (PentaRay, Biosense Webster) having a five-branched star design which allows for high-density mapping, it is shown in Figure 2.8. Non-contact Mapping: The EnSite Velocity System (St. Jude Medical Inc.), shown in Figure 2.9, is part of the non-contact mapping systems. In this approach, the recording device consists of a catheter with a multi-electrode array of 64 unipolar electrodes over an inflatable balloon, each of which can remotely measure the potential generated by far-field electrograms [Aliot 2009]. Because the electrical activity does not correspond to the potentials at the endocardial surface but in the blood pool, an inverse problem must be solved. Virtual unipolar electrograms at the endocardium are accurately reconstructed and in agreement with contact mapping recordings if the surface of interest is within 4 cm of the center of the MEA balloon [Gornick 1999] [Schilling 1998] [ Sivagangabalan 2008]. Because of these characteristics, the system is preferred for activation mapping. It also has the advantage of being applicable in patients who cannot tolerate arrhythmia or when an arrhythmia cannot be reproduced during the EP study.

Electrogram Interpretation

Aside from the mapping system chosen for a given procedure, it is of vital importance for the clinical electrophysiologist to be able to recognize complex patterns in recordings inherent to the nature of the electrograms (EGM), in order to accurately diagnose and locate arrhythmias [Stevenson 2005].

Cardiac electrograms are the voltage differences between two recording electrodes (a positive lead, known as the anode, and a negative lead, referred to as the cathode) at a given point in the cardiac cycle [Stevenson 2005]. There exist two main kinds of electrograms recorded during EAM: unipolar and bipolar.

Unipolar electrograms: When recording unipolar electrograms, the anode is in contact with the myocardial tissue of interest and the cathode is placed at a distant point from the heart, such that its influence in the cardiac measurements is negligible. Although theoretically the cathode should be placed at an infinite Chapter 2. Clinical and Technical Context distance from the heart to obtain a truly unipolar recording, in practice, it is located in the inferior vena cava [Stevenson 2005]. Following this convention of electrode placements, a depolarization wavefront moving towards the electrode in the cardiac chamber will result in a positive deflection in the EGM, whereas a receding one will yield a steep negative deflection. The point of maximum negative slope (-dV /dt) corresponds to a wavefront directly underneath the recording electrode [Stevenson 2005]. This characteristic renders unipolar electrograms valuable for the generation of activation maps and to identify the site of earliest activation in focal tachycardias. One major limitation of unipolar electrograms is the effect of far-field signals from depolarising myocardial tissue at the distance [Aliot 2009], this is particularly hindering in scar regions where the local low-voltage potentials can be overshadowed by far-field effects. Sample unipolar electrograms are shown in the top and middle rows of

Electro-anatomical mapping

are generated by computing the difference in the potential recorded between two closely-spaced electrodes. Local activation is identified by the initial signal peak, but this measure is less reliable than with unipolar electrograms if the signals are fractionated or have multiple components. The far-field signals in the two electrodes is considered equivalent, therefore bipolar EGMs are useful to evaluate local signals [Stevenson 2005], even those with low-amplitude potentials, making them ideal for mapping scar-related arrhythmias [Aliot 2009].

Most catheters are equipped with four electrodes numbered from 1 to 4 from the distal-most to the proximal-most. Comparisons between the bipolar signals generated by neighboring electrodes indicate wavefront direction, but the amplitude of the signals can be severely compromised if the wavefront is moving perpendicular to the axis of the recording dipole.

In order to better exploit the strengths of both types of electrograms, EP studies usually obtain both unipolar and bipolar recordings simultaneously.

Mapping Approaches

Different mapping approaches can be used to study a particular type of arrhythmia in a given patient, though they are often complementary. This section describes the most commonly used ones:

Activation Mapping: Consists in the recording of the electrical activation sequence of the chamber of interest. It is necessary to establish both a location and timing reference, as well as a timing window. Although in theory either intracardiac EGMs or surface ECG leads can be used to chose a timing reference, intracardiac EGMs are preferred [START_REF] Bhakta | [END_REF]].

Pace Mapping: During this approach, a site of interest is paced during the absence of VT and the resulting activation sequence and surface ECG morphology are assessed. If the pacing of a particular site reproduces the same ECG morphology as the one seen during clinical VT, then the location is suspected to be either the origin of a focal VT or the exit of a scar-related re-entry [START_REF] Brunckhorst | [END_REF]].

Entrainment Mapping: Is performed by delivering a train of stimuli at a rate just higher than the one present during VT [Stevenson 1993]. It serves to identify re-entry circuits sites, e.g. surface ECG morphology during entrainment will determine if the pacing site is within an isthmus or not [Aliot 2009].

Substrate Mapping: Refers to the characterization of a region of interest, particularly those suspected to promote re-entry due to their anatomical or electrophysiological characteristics. Scarred tissue is frequently the target of substrate mapping and its correlation with low-voltage areas has been found in animal models [Callans 1999]. Regions with particularly low-voltage (<0.5mV) have been referred to as dense scar but can still participate in re-entry circuits as they can contain varying amounts of viable myocytes [Soejima 2002]. Furthermore, targeting the border zone for ablation, which consists of regions with voltage ranging from 0.5 to 1.5 mV, has gained interest in the recent years [Marchlinski 2000] [Aliot 2009].

Radiofrequency Catheter Ablation

Radiofrequency ablation (RFA) is a technique in which a tissue region is subjected to extreme heat resulting from a high-frequency alternating current [Townsend Jr 2012]. In cardiology, it is used to destroy electrical pathways responsible for cardiac arrhythmias. The extent of the RFA lesion created is dependent of operator-controlled characteristics, such as RF power and duration (which should be at least 35 to 45 seconds to approach steady state [Haines 1993]), but also on operator-independent phenomena, including electrode-tissue contact and the cooling effect by the circulating blood [Haines 1993]. Temperature monitoring is very important in order to ensure a proper ablation of the target site while minimizing the risk of serious complications such as of thromboembolism, although these are infrequent [Matsudaira 2003].

RFA is considered the sole treatment option for patients suffering from recurrent ventricular tachycardias (VT) without structural heart disease [Aliot 2009]. Patients with scar-related VTs are also referred for RFA therapy in conjunction with implantable-cardioverter defibrillators (ICDs) and/or anti-arrhythmic drug therapy [Aliot 2009]. RFA ablation is a complex intervention reserved to patients with advanced heart disease [Aliot 2009]. The percentage of significant complications (requiring prolonged hospitalization, additional interventions or death) reported in [Reddy 2007] was of < 5% on patients undergoing prophylactic catheter ablation of post-infarction VT in a randomized, multi-center study.

Ablation Target Identification Up Until Now

Though RFA has proven to be an effective therapy in patients suffering from ventricular tachycardia (VT) whose ICD triggers frequent defibrillation shocks [Aliot 2009], to date, there exists no universal consensus on the optimal ablation strategy.

A number of techniques are found in the literature to treat either unmappable or hemodynamically intolerable VT. Studies such as [Marchlinski 2000] reported the use of linear endocardial RFA lesions going from the regions of dense scar to those of normal myocardium to control unmappable VT. [Kottkamp 2003] described how 80% of the 28 patients in their study were rendered non-inducible after the placement of line lesions that transected all potential isthmuses. The study in [Soejima 2001] also used lines of RF lesions to treat unstable VT and suggested that less ablation is needed when a reentry circuit isthmus is identified.

Furthermore, reentry isthmuses can be identified in a number of ways. In [Soejima 2002], these were located by identifying electrically unexcitable scar within low-voltage infarct regions. Other strategies include the assessment of pace maps,

Radiofrequency Catheter Ablation

where an abrupt change in paced QRS morphology is hypothesized to be an index of a successful ablation target in patients with well-tolerated post-infarct endocardial re-entrant VT [De Chillou 2014] or through pace-mapping within the scar tissue with the aim of identifying electrogram characteristics that are helpful for isthmus identification during sinus rhythm, as done in [START_REF] Bogun | [END_REF]]. The authors explained how good pace maps at the sites of isolated potentials within the scar identify parts of the critical isthmuses in post-infarction VT patients.

Voltage maps and local electrogram voltage characteristics have also been studied in the identification of slow conduction channels responsible for re-entrant circuits. The study in [Hsia 2006] defined a conducting channel as a path of multiple activated sites within the VT circuit that demonstrated an electrogram amplitude higher than that of surrounding areas and concluded that most entrance and isthmus sites of hemodynamically stable VT are located in dense scar, whereas exits are located in the border zone. Nonetheless, these approaches require a careful voltage threshold adjustment. In [Arenal 2004], the effect of different levels of voltage scar definition in the identification of conducting channels inside the scar was studied in a patient cohort with chronic myocardial infarction referred for VT ablation. Although most of the conducting channels were identified by defining the voltage scar with a threshold of ≤ 0.2mV and RFA supressed the inducibility in 88% of the tachycardias, they concluded that a tiered decreasing-voltage definition of the scar is critical for channel identification. [START_REF] Biase | [END_REF] suggested an aggressive scar homogenization technique in which all abnormal electrograms were targeted, with successful results in endo-epicardial applications. Others, like [Tilz 2014] described how the use of linear ablations encircling the scar would yield its electrical isolation and had promising results in VT treatment.

Endpoints for RFA Therapy

Traditionally, VT non-inducibility is still considered the endpoint in RFA treatment. According to [Aliot 2009], three general endpoints to catheter ablation have been evaluated:

1. Non-inducibility of clinical VT 2. Modification of induced VT cycle length (elimination of all VTs with cycle lengths ≥ spontaneously documented or targeted VT)

3. Non-inducibility of any VT (excluding ventricular flutter and fibrillation)

Nonetheless, these are not guaranteed to be the optimal RF ablation endpoints and debate still exists in their use for assessing ablation lesion creation [Aliot 2009].

Alternative Endpoints to RFA Therapy: In recent years, new endpoint strategies have been proposed. For the treatment of VT in 18 patients with arrhythmogenic right ventricular cardiomyopathy (ARVC), [Nogami 2008] evaluated the modification of RFA therapy endpoint to the change in isolated delayed components scar isolation [Tilz 2014] Non-exhaustive review of proposed strategies for substrate-guided VT ablation. Modified from [Berruezo 2014] (IDC). They concluded both that IDCs that appear in sinus rhythm are related to clinical VT and can therefore be used as ablation targets and that a change in IDC was a possible endpoint to the therapy due to the siginficantly lower VT recurrence in patients who presented it with respect to those with no IDC or unchanged IDC.

A recent study with a larger cohort of 70 patients suffering from VT and structurally abnormal ventricles proposes the elimination of local abnormal ventricular activities (LAVA) as an endpoint for RFA therapy [START_REF] Jaïs | [END_REF]]. The clinicians ablated LAVA found through the use of a high-density mapping catheter during both sinus rhythm or ventricular pacing and concluded that their elimination was associated with a reduction in recurrent VT or death during long-term follow-up [START_REF] Jaïs | [END_REF]]. Table 2.1 presents a non-exhaustive review of proposed strategies for substrateguided VT ablation.

Cardiac Electrophysiology Simulation

There exist several models for computing cardiac electrophysiology ranging from highly complex to fairly simple. They can be broadly classified into three groups: biophysical models, phenomenological models and Eikonal models.

-Biophysical models deal with a highly complex modelling of the different ionic concentrations and channel dynamics at a cellular level and can include over 50 variables; an example of them is the ten Tusscher model [K. ten Tusscher 2004].

-Phenomenological models often appear as simplifications of the biophysical models, with a consequent reduction in the number of involved variables. Nonetheless, they are still able to capture the shape of the action potential and are useful to model the propagation at the organ scale. The Aliev-Panfilov [Aliev 1996], Fenton-Karma [Fenton 1998], FitzHugh [FitzHugh 1961] and [START_REF] Mitchell | A Two-Current Model for the Dynamics of Cardiac Membrane[END_REF]] models belong to this group. These models can be further classified into bi-domain or mono-domain. Bi-domain models consider an intracellular and an extracellular domain separated by the cell membrane but which remain overlapping and continuous. It models the evolution of the potentials in both domains separately. Mono-domain models, on the contrary, assume the extracellular potential is grounded. Therefore, the intracellular potential is equal to that of the membrane.

-Eikonal models [Keener 1991] correspond to static non-linear partial differential equations of the depolarisation time derived from the previously mentioned models. They are not able to account for complex physiological states.

These models have been used for numerous applications [Trayanova 2011] varying from the study of the normal wave propagation in the ventricles [START_REF] Simelius | [END_REF]] and the effect of repolarization parameter heterogeneity across the heart [Franzone 2008] to the understanding of the organization of human ventricular fibrillation to identify potential therapeutic targets [START_REF] Tusscher | Organization of ventricular fibrillation in the human heart: experiments and models[END_REF].

The Mitchell-Schaeffer Model

Throughout this thesis, the Mitchell-Schaeffer (MS) model was used to model cardiac electrophysiology. It is a simplified biophysical model derived from the Fenton-Karma ionic model [Mitchell 2003] [Fenton 1998]. It has found applications in patient-specific personalisation for VT simulation [Relan 2011a] and interactive simulator generation of patient-specific electrophysiology [Talbot 2014].

It models the transmembrane potential as the sum of a passive diffusive current and several active reactive currents including a combination of all inward and outward phenomenological ionic currents. This is described in the MS model by a system of partial differential equations:

       ∂ t u = div(D∇u) + zu 2 (1-u) τ in -u τout + J stim (t) ∂ t z = (1-z) τopen if u < u gate -z τ close if u > u gate (2.1)
Where:

u is a normalized transmembrane potential variable z is a gating variable which depicts the depolarization and repolarization phases by opening and closing the currents gate

AP D n+1 = f (DI n ) = τ close ln( 1 -(1 -h min )e -DIn τopen h min ) (2.2)
where h min = 4( τ in / τout ) and n is the cycle number. The maximum value of APD is also explicitly derived as:

AP D max = τ close ln( 1 h min ) (2.3)
CV restitution. This formulation is also present in the MS model. Its mathematical formulation is described in the following equation, with g(DI) = CV as the next cycle CV:

g(DI) = ( 1 4 (1 + 1 - h min h(DI) ) - 1 2 (1 -1 - h min h(DI) ) 2dh(DI) τ in ) (2.4)

Machine Learning in Medical Imaging

Deriving realistic mathematical models to explain complex physiological phenomena can sometimes prove to be a difficult task. Machine learning is a branch of artificial intelligence whose aim is the generation of algorithms capable of learning from observations in order to make predictions about future ones. It has found applications in fields as varied as language processing [Collobert 2008] [Daelemans 2003], computer vision [START_REF] Taigman | [END_REF]] [Huang 2008] or stock market analysis [Huang 2005] [Schumaker 2009]. Its increased use by the medical imaging community is due to the processing of high volumes of multi-dimensional data, including uncertainty and noise, coming from multi-modal sources that exist in this domain. Nonetheless, one of the main challenges is the generation of ground truth, as it requires the annotation of large volumes of data by a set of medical experts. Therefore, large efforts are continuously being made to provide clinicians with user-friendly and efficient annotation tools. Among the various applications of machine learning techniques in the medical field are detection or segmentation of anatomical structures [Criminisi 2011a] [Geremia 2013], registration [Muenzing 2012] [Aljabar 2012] and image recognition [START_REF] Margeta | [END_REF]].

Classification tasks are among the functions that machine learning techniques are used for and many of such algorithms exist in the literature. One of the most widely used is the support vector machines (SVM) classifier because it guarantees maximum-margin separation in binary classification tasks. Boosting also ranks among the most popular techniques and works by building a strong classifier as a linear combination of many weak ones. Random forests have emerged as a machine learning technique that can simply and effectively combine randomly trained classification trees and has been shown to be easily extendable to multiple class problems [Criminisi 2011a].

Random Forests for Classification

Random forests are discriminative classifiers that generate a hierarchical representation of the training data which is optimized for testing. They have been chosen as the machine algorithm framework of this thesis because, unlike state-of-the-art classifiers such as SVMs, they are created in an intuitive and easily understandable structure. Moreover, they also provide informative uncertainty measures on the classification results [Criminisi 2011a], which would only be possible with probabilistic SVM implementations such as the relevance vector machine at the expense of computational efficiency. They also possess maximum-margin properties, the hallmark of SVMs, implicitly incorporate the notion of bagging into their framework and compare to boosting techniques by building strong classifiers from a large number of weaker ones. This framework has already successfully found multiple applications in the medical image processing community, such as brain lesion segmentation [Geremia 2013], myocardium delineation in 3D echocardiography [START_REF] Lempitsky | Random forest classification for automatic delineation of myocardium in real-time 3D echocardiography[END_REF]], segmentation of the left atrium in 3D MRI [Margeta 2014b] and organ localization in CT volumes [Criminisi 2011a].

The goal in classification is to automatically associate an input data point v with a discrete class c ∈ {c k }. Machine learning algorithms have proven to be powerful tools to discover the relevant correlations within a dataset and to design a classification method based on the most discriminant features. The training data consists of a set of labeled observations T = v k , Y (v k ) where the label Y (v k ) is most commonly given by an expert clinician in the medical imaging domain. When asked to classify a new input observation, the classifier aims to assign it a label y(v). The process of classification can be divided into an offline training phase and an online testing phase, both of which will be described next.

Training phase

During training, exemplified in Figure 2.11, the algorithm optimizes the parameters of the weak learner at each split node (j):

θ * j = argmax I j θ j ∈τ j (2.5)
where the objective function I j is the information gain for discrete distributions:

I j = H(S j ) -Σ i∈L,R | S i j | | S j | H(S i j ) (2.6)
with i ∈ L, R representing the left and right child nodes. The term H(S) corresponds to the Shannon entropy:

H(S) = -Σ c∈C p(c) log p(c) (2.7)
and p(c) is calculated as the normalized empirical histogram of labels corresponding to the training points in S. It can be seen in Figure 2.11 that the maximization of the information gain in the trees produces entropy of the class distributions that decrease, and therefore the prediction confidence increases, as you move towards the leaves.

A degree of randomness is injected to the trees during this phase by random training data set sampling (bagging), which corresponds to providing the tree with only a partition of the training data as input [Criminisi 2011a]. It is important to note that the generation of the random subset τ j can be done on the fly.

The choice of parameters in the generation of the forest has an impact in the model predictions, particularly that of the forest size and the tree depth [Criminisi 2011a]. Their effects can be briefly summarized as following:

-Forest Size. The forest has T components with t indexing each tree. Single trees produce over-confident predictions, which are likely to yield erroneous classifications when presented to test data. With an increase of the number of trees and the averaging of their posteriors, forests improve their generalization behavior by producing smoother posteriors, as showin in Figure 2.12.

-Tree Depth. This parameter controls the amount of overfitting. Forests built with trees that are too shallow produce low-confidence posteriors, whereas increasing the depth of the forests can lead to overfitting.

Testing phase

Once the trees are fixed during the training phase, the testing phase (depicted in Figure 2.13) is completely deterministic. When applied to a new test signal

p(y( v) = b) = 1 T T t=1 p lt (Y ( v) = b) (2.8)
It is important to highlight that, as shown by the previous equation, classification forests produce a probabilistic output by returning not only a single class point prediction but an entire class distribution [Criminisi 2011a].

Conclusion

In this chapter, we presented a background of cardiac anatomy and electrophysiology, in both the healthy heart and the pathological arrhythmia scenario. We also described the technologies being used to date in clinics to map cardiac electrophysiology in a minimally invasive manner. Radio frequency ablation as a potentially curative treatment was presented as well as the challenges it presents. On the technical side, a description of the types of mathematical models being used in cardiac eletrophysiology research was included. Lastly, we briefly introduced the growing field of machine learning in the medical domain.

Introduction

Current risk stratification in patients who are at risk of potentially fatal ventricular arrhythmias but without a prior history of sustained arrhythmia relies on determination of left ventricular (LV) function, the presence of myocardial scar and arrhythmia inducibility during electrophysiological testing [Moss 1996[START_REF]TA095: arrhythmiaâimplantable cardioverter defibrillators (ICDs) (review)[END_REF]]. In patients with ventricular arrhythmias, especially those with ischaemic cardiomyopathy (ICM), radiofrequency ablation (RFA) is increasingly used to treat ventricular tachycardia (VT) in order to reduce implantable cardioverter defibrillator (ICD) discharges, improve patient quality of life and reduce mortality as ICD shocks are a cause of substantial morbidity. The current risk stratification strategy is imperfect with not all high-risk patients receiving an ICD and those receiving one never experiencing appropriate therapies. Similarly, ablation of VT is technically challenging with a recurrence rate of up to 40% with a lack of clinical consensus on the optimal ablation strategy [Aliot 2009]. Better risk stratification and higher ablation success rates would potentially improve patient outcomes. There is therefore a need to identify individuals at high risk of developing ventricular arrhythmia and the arrhythmia substrate amenable to RFA in order to guide the optimal ablation strategy [Stevenson 1993].

Computational modelling of cardiac arrhythmogenesis and arrhythmia maintenance have made a significant contribution to the understanding of the underlying mechanisms of arrhythmia [START_REF] Courtemanche | [END_REF]] [Watanabe 2001[START_REF] Panfilov | [END_REF]] [Jalife 1996] [Cherry 2004] [Clancy 1999]. Studies have identified multiple factors involved in the onset of arrhythmia, including wave fragmentation, spiral wave breakups, heterogeneity of repolarization, action potential duration (APD) restitution and conduction velocity (CV) [START_REF] Trayanova | [END_REF][START_REF] Killeen | [END_REF]] [Pueyo 2011] [Yue 2005] [Arevalo 2007[START_REF] Banville | [END_REF]]. In particular, the heterogeneity in APD restitution, the adaptation of APD as a function of the cardiac cycle length, has a crucial role in arrhythmogenesis [Cherry 2004[START_REF] Nash | Whole heart action potential duration restitution properties in cardiac patients: a combined clinical and modelling study[END_REF]] [Clayton 2011]. Imagebased computational models have incorporated cardiac structural information into such simulations [Relan 2011a] [Ashikaga 2013]. However the integration of both personalized structural and functional data has not previously been performed.

We hypothesized that using electrophysiological mapping data and structural anatomical data acquired respectively from invasive electrophysiology studies and high-resolution cardiac magnetic resonance imaging (MRI), we could develop a patient-specific biophysical model to evaluate how these properties were involved in the induction of VT. In a subgroup of the patients, we also compared in silico VT stimulation studies against clinical VT stimulation studies conducted in the cardiac 3.2. Methods 29 catheterization laboratory to predict both VT inducibility and the characteristics of the induced VT circuits.

Methods

Patients

All patients were prospectively invited to participate in the study following local research ethics committee approval and all patients gave written consent prior to study inclusion. We studied 7 patients, 5 with ischaemic cardiomyopathy (ICM) and 2 with non-ischaemic dilated cardiomyopathy (NICM) defined in accordance with the criteria of the World Health Organisation/International Society and Federation of Cardiology [Richardson 1996]. All patients were being considered for ICD implant for primary prevention on the basis of their LV function. All patients underwent MRI with delayed enhancement (DE) for scar assessment prior to their electrophysiology study to acquire cardiac anatomical and functional information. Three patients underwent a clinically indicated VT stimulation study.

MRI Acquisition and Image Processing

All 7 patients completed MRI morphological and volumetric assessment as well as scar characterization by DE-MRI. Imaging was performed on a Philips Achieva 1.5T scanner using a 32 channel cardiac coil. A high-resolution three-dimensional (3D) whole heart balanced steady-state free precession (B-SSFP) free-breathing scan with (acquired) isotropic resolution of 1.8 mm 3 was performed using respiratory navigator motion correction for the purpose of cardiac structure segmentation. LV volumetric assessment was performed using a standard stack of short-axis B-SSFP slices. High-resolution scar imaging was acquired using a free-breathing respiratory navigated inversion-recovery sequence 20 minutes post intravenous injection of a gadolinium contrast agent (Gadobutrol 0.2mmol/kg), with an acquired voxel size 1.3×1.3×2.6mm 3 , field of view approximately 300×300×100mm, repetition time/echo time of 5.4/2.6ms and flip angle 25 degrees. Patient-specific inversion time for the sequences was selected individually based on a preceding Look-Locker scan to ensure the optimal nulling of the myocardium.

LV volume and ejection fraction were calculated from the short-axis B-SSFP images using a ViewForum Workstation (Philips Healthcare, Netherlands). The 3D B-SSFP structure images were processed to obtain a structural model of the four chambers of the heart using software developed within an open source frame-work GIMIAS [Peters 2007[START_REF] Larrabide | [END_REF]]. The LV myocardial scar distribution was segmented using signal intensity (SI) based analysis from the high-resolution DE-MR images. Using the full width-half-maximum (FWHM) method, all voxels with SI values above the half of the maximum SI were automatically characterized as scar core. The standard deviation of a manually selected remote region of presumed non-infarct myocardium was computed. Pixels with SI higher than twice this stan-Chapter 3. VT Inducibility Prediction: A combined modeling and clinical approach dard deviation (2SD) but lower than that identified as scar core were automatically assigned as gray zone (an admixture of scar and healthy myocardium scar, often in the region of scar border zone) [START_REF] Kim | [END_REF]]. Finally, a personalized 3D model of the ventricles was derived from the MRI images: a tetrahedral mesh was generated from the binary mask of the ventricles using the CGAL software (http://www.cgal.org). Each element of the mesh was labelled (healthy / scar core / gray zone) according to the segmentation of the myocardium performed in the previous step.

Electroanatomical Mapping and Signal Processing

LV non-contact electroanatomic mapping (EAM) was performed using a multielectrode array catheter (EnSite Velocity System, St Jude Medical, MN, USA) passed via the femoral artery retrogradely across the aortic valve into the LV cavity in all seven patients. The EnSite system uses the inverse solution method to reconstruct endocardial electrical potentials within the LV cavity [Schilling 1998].

The chamber geometry was reconstructed using locator signals from a steerable electrophysiological catheter. Three patients (Patients 1-3) with ICM underwent a simultaneous VT stimulation study according to the Wellens protocol with pacing from the RV apex during the mapping study [Wellens 1985] as part of the clinical work-up for risk stratification to determine if an ICD should be implanted as per National guidelines [START_REF]TA095: arrhythmiaâimplantable cardioverter defibrillators (ICDs) (review)[END_REF]]. Unipolar electrograms (UEG) derived from non-contact mapping were filtered from an electrophysiology recorder (EnSite Velocity System) with a band-pass filter. In order to optimize QRS complex and T wave detections, the high-pass and low-pass filter cut-off frequencies were set respectively at 10Hz/300Hz and 0.5Hz/30Hz. The data were then exported for offline analysis. The depolarization times were detected within the QRS window and derived from the zero crossings of the laplacian of the measured UEGs [START_REF] Coronel | [END_REF]]. The repolarization times were detected within the ST window for the signals and derived using the alternative method [START_REF] Arthur M Yue | Determination of Human Ventricular Repolarization by Noncontact Mapping Validation With Monophasic Action Potential Recordings[END_REF]]. The alternative method has repolarization times derived from -dV /dt max for the negative T-wave, at the dV /dt min for the positive T-wave, and the mean time between -dV /dt max and dV /dt min for the biphasic T-waves.

The relationship between the diastolic interval (DI) of one cardiac cycle and the APD of the subsequent cardiac cycle is described by the APD restitution curve (APD-RC). The difference between the depolarization time and repolarization time was used to estimate the activation recovery time (ARI), which is a surrogate marker for APD. The APD-RC was estimated during steady state RV pacing (600ms, 500ms and 400ms) with sensed extras at different coupling intervals. The APD-RC was represented by a non-linear equation using a least-squares fit to the mono-exponential function as previously detailed on experimental and clinical data [Relan 2011a] [Relan 2011b]: a single APD-RC was fitted for each measured point from the EAM and the maximum APD was estimated as the asymptotic APD of the APD-RC, when the DI tends to infinity.

Methods

Cardiac Electrophysiology Models

The present study used a coupled personalization framework that was previously evaluated in details [Relan 2011a] [Relan 2011b]. It combines the benefits of two different kinds of mathematical models while keeping the computational complexity tractable. The Eikonal (EK) model was used to estimate the conductivity parameters over the ventricle derived from non-contact mapping of the ventricular endocardial surface potential, which were then used to set the parameters for the Mitchell-Schaeffer (MS) model. Additionally, the MS model is able to hold the memory of one preceding cycle and has restitution properties, thus is able to simulate arrhythmias macroscopically. Both models are computed on the whole 3D myocardium. This personalization framework has already been detailed in a previous publication and the predictive power of such personalized model was evaluated on experimental data [Relan 2011a] [Relan 2011b]. The process of building the models from the MRI and EAM data is illustrated in Figure 3.1 and a summary on the used models and their personalization are included in the following sections.

c o d EK ( ∇T (x) t M ∇T (x) -∇ • (d EK M ∇T (x))) = τ (x) (3.1)
Where :

c o is a dimensionless constant (=2.5) τ (x) is the cell membrane time constant (=0.003s)
d EK is the square of the tissue space constant along the fiber and is related to the specific conductivity of the tissue in the fiber direction, it has units of m 2 M = diag(1, ρ, ρ) is a diffusion tensor that incorporates the tissue anisotropy into the model ρ = 1 / 2.5 2 is the anisotropy ratio between longitudinal and transverse diffusion. This value is given to have a conduction velocity (CV) 2.5 times faster in the fiber direction Details for solving the previous non-linear equation are explained in a prior publication [START_REF] Sermesant | [END_REF]].

Mitchell Schaeffer (MS) Model

The MS model is a simplified biophysical model derived from the Fenton Karma ionic model [Mitchell 2003] [Fenton 1998]. It models the transmembrane potential as the sum of a passive diffusive current and several active reactive currents including combination of all inward and outward phenomenological ionic currents. This is described in the MS model by a system of partial differential equations:

       ∂ t u = div(D∇u) + zu 2 (1-u) τ in -u τout + J stim (t) ∂ t z = (1-z) τopen if u < u gate -z τ close if u > u gate (3.2)
Where:

u is a normalized transmembrane potential variable z is a gating variable which depicts the depolarization and repolarization phases by opening and closing the currents gate

J in = zu 2 (1-u)
τ in represents the inward currents (primarily N a + and Ca 2 +) which raise the action potential voltage J out = u τout represents the outward currents that decrease the action potential voltage (mainly K + ), describing repolarization

Methods

33

J stim is the stimulation current at the pacing location τ in , τ out , τ open , τ close have units of seconds M is a diffusion tensor that controls the diffusion term in the model This model incorporates both APD and CV restitution effects, and the restitution curves can be written in an analytical formulation, which can be used to estimate restitution parameters.

Model Personalization

Apparent Conductivity Parameter Estimation. The apparent conductivity (AC) of tissue can be measured by the parameter d EK in the EK model. It is initially estimated on the endocardial surface as a global value using a bisection method that matches the average conduction velocity of the measured depolarization time (DT) isochrones to the simulated ones. It is used as an initial estimate in an adaptive multi-level domain decomposition algorithm, which minimizes the mean-squared difference of the simulated and measured DT isochrones [Chinchapatnam 2008]. One particular aim of this step is to estimate areas with high conduction velocity, which can represent parts of the Purkinje network.

Extrapolation of myocardial volumetric parameters from endocardial data, imaging and body surface ECG In order to propagate the estimated parameters to the whole myocardium, the present approach uses a diffusion-like process that smoothly extrapolates the estimated endocardial parameters. This estimation is based on the model using the anatomical 3D information from the MRI scar images and the global electrical information (ECG). The transmural AC value is estimated in order to fit the measured QRS duration, by a one-dimensional minimization of the following cost function: (mean-squared difference of simulated and measured isochrones at endocardium + squared difference of simulated and measured QRS duration). The simulated QRS duration is calculated as the difference between the maximum and the minimum depolarization times in the biventricular mesh and the measured QRS duration is estimated from the surface ECG. The AC values for RV endocardium and RV myocardial mass are set at 5 mm 2 and 0.64 mm 2 as previously described [START_REF] Keener | Mathematical physiology: I: Cellular physiology[END_REF]]. Such estimation method has been validated by experimental data [Relan 2011b]. Due to the absence of transmural electrical propagation information, no variation is assumed across the healthy left ventricular myocardium, excluding the LV endocardium which was personalized as described in the previous paragraphs and the scar personalized from imaging, see below.

Coupling of EK and MS Model Parameters

The AC parameter for d EK model is a scale for the diffusion speed of the depolarization wave front in the tissue. In Chapter 3. VT Inducibility Prediction: A combined modeling and clinical approach 3D, the model Conduction Velocity (CV) is related to d EK as :

c EK = α EK d EK + β EK (3.3)
where the constants α EK and β EK are introduced to take into account the discretization errors (in particular of the curvature) in 3D.

The corresponding conductivity parameter for MS model, d M S is also a scale for the wave diffusion speed in the tissue. The model CV here is related to d M S as:

c M S = α M S d M S + β M S (3.4)
where the constants α M S and β M S are introduced for the same reasons as of EK model, while τ in is kept as a constant. The estimated AC parameter d EK can then be used to estimate the parameter d M S as previously described [Relan 2011a].

To represent the volumetric heterogeneity of the border zone, the intensity values derived from the 3D DE-MRI images were used to modulate the spatial variation in the AC values of scar border zone (the higher the signal intensity seen on DE-MRI images, the lower the tissue conductivity). AC values were assigned to zero in the scar core, to reflect the paucity of viable cardiac myocytes.

APD Restitution Curves. APD restitution is an electrophysiological property of the cardiac tissue and defines the adaptation of APD as a function of the heart rate. Its slope has a heterogeneous spatial distribution. The APD restitution curve (APD-RC) defines the relationship between the diastolic interval (DI) of one cardiac cycle and the APD of the subsequent cardiac cycle. The slope of these RCs is controlled by a model parameter τ open of the MS model and depicts the APD heterogeneity present at multiple heart rates. APD-RC for MS model is explicitly derived as:

AP D n+1 = f (DI n ) = τ close ln( 1 -(1 -h min )e -DIn τopen h min ) (3.5)
where h min = 4(τ in /τ out ) and n is the cycle number. The maximum value of APD is also explicitly derived as AP D max = τ close ln(1/h min ).

CV Restitution Curves. In the current study, a greater emphasis is placed on the impact of APD restitution, as in our experience CV restitution has a greater impact on propagation at a much higher pacing frequencies than our sampling pacing range. If deemed relevant, CV restitution can also be personalized, as the parameters to adjust CV restitution are also present in the Mitchell-Schaeffer model. Its mathematical formulation is described in the following equation: with g(DI) = CV as the next cycle CV.

g(DI) = ( 1 4 (1 + 1 - h min h(DI) ) - 1 2 (1 -1 - h min h(DI) ) 2dh(DI) τ in ) (3.6)

Methods

However, these parameters also impact the APD restitution formula, so a slightly more complex model like the Fenton-Karma model may be better suited in such application [Fenton 1998].

In silico VT Stimulation Study

In 3 patients with ICM, clinical VT stimulation studies were carried out at the same time as EAM. In these patients, the personalized model was used to simulate a clinical VT stimulation study in silico in order to predict the initiation and maintenance of ventricular arrhythmia. The simulated study was performed in accordance with the clinical study protocols that were carried out in these patients: the pacing stimuli were applied from the RV apex following the Wellens protocol. If a sustained VT was initiated in silico, the VT cycle isochrones were computed to locate the site of the earliest activation corresponding to the exit point of the re-entrant VT circuit. These results were compared with the observed results from the clinical VT stimulation studies. Each stage for the VT stimulation in silico was simulated in parallel on a cluster of computers. The personalization plus VT simulation study for a mesh of >230,000 tetrahedra runs in approximately 12 hours.

In order to test in silico the inducibility of VT from other pacing sites, the VT stimulation study was also simulated from alternate pacing sites in the RV and LV that encompassed the basal and apical freewall/lateral/septum/anterior and inferior walls. This allowed identification of all the potential VT circuits and exit points of induced sustained monomorphic VT from the simulated study. The exit points identified during the clinical VT stimulation studies as well as the exit points identified during the in silico simulations were characterized in terms of the spatial heterogeneities of the AC, APD restitution properties, scar and gray zone at their locations.

VT simulation with solely image-based personalized parameters. A sensitivity analysis where the performance of a model constructed using non-personalized empirical electrical parameters was compared to that using personalized patientspecific electrical data. For all cases, patient anatomy and scar geometry was preserved. The VT simulation procedure was performed as described in the previous sections.

The personalized mesh was divided into three different types of tissue, guided by the delayed enhanced MR imaging information: healthy myocardium, grey zone and scar core. Nodes within the scar core tissue were set to have zero apparent conductivity values. Conduction and restitution parameters in the healthy myocardium tissue were set to the average value computed from the personalized models across patients. This yielded values of 4 m 2 s 2 for the apparent conductivity and τ close = 180ms, τ in = 0.3ms, τ out = 6ms, τ open = 160 ms for the restitution model parameters.

According to the described values of conduction velocity and restitution parameters in the literature, model parameters were modified in the grey zone areas [Yao 2003] [Decker 2010]. Conduction velocity parameters were reduced by 90% clinical approach with respect to those in the healthy myocardium. Restitution parameter values in the grey zone regions were set to have a 25% longer action potential duration (increase in τ close parameter) and a 31% smaller upstroke velocity (increase in the τ in parameter) than those in the healthy myocardium.

Statistical Analysis

The electrical parameters consisted of multiple data points, each corresponding to a different location within the LV of the individual patients. The data were expressed by median and inter-quartile range (IQR) or mean and standard deviation. Continuous variables were compared using Median test. A p value < 0.05 was considered to be statistically significant. All statistics were performed using computer software SPSS Statistics, version 21, SPSS Inc., IBM, USA.

Results

Table 3.1 summarizes patient characteristics and MRI findings. During the clinical VT stimulation studies, patients 1 and 2 developed sustained monomorphic VT, while patient 3 was non-inducible. All seven patients were implanted with an ICD for primary prevention based on guidelines following clinical assessment. During a median follow-up period of 22 months, IQR 9 months, Patient 2 and Patient 4 received appropriate ICD therapy for sustained VT.

Spatial Heterogeneities in AC and APD Restitution Properties

The spatial heterogeneities of the LV myocardium in terms of apparent conductivity and APD restitution properties were characterized for each patient. The estimated APD restitution properties and tissue conductivity across the LV endocardium are depicted in Figure 3.2. The range of AC across the LV was comparable across the 7 patients. The median APD-RC slopes and maximum asymptotic APD for NICM patients was lower than that for ICM patients (p<0.01). Amongst the three patients who underwent a clinical VT stimulation study, Patients 1 and 2 with positive clinical VT stimulation studies had APD-RC slopes IQRs of 1.50 and 1.02 and maximum APD IQRs of 123ms and 66ms respectively. The range of APD-RC slopes and maximum asymptotic APDs were greater than that of Patient 3 with a negative clinical VT stimulation study who had an APD RC slope IQR of 0.48 and a maximum APD IQR of 36ms. Patient 4 with appropriate ICD therapy on follow-up had a greater APD RC slope IQR than Patients 5-7 with no appropriate ICD therapy on follow-up.

The APD-RC properties and ACs are illustrated on LV polar plots in Figure 3.3 to provide a qualitative assessment of the spatial heterogeneity of these parameters across different regions of the LV. There appeared to be a broad spatial heterogeneity in AC in all the study patients including the two patients with NICM who also had regions of low AC in areas without DE seen on MR scar images. However, there was a wide variation in the spatial heterogeneity of APD-RC properties across the cohort, with patients 1 and 2 showing the most heterogeneity. The VT exit points that were observed during the clinical VT stimulation studies in Patients 1 and 2 appear to be at the region with the steepest gradients in heterogeneity of APD-RC properties and AC. In addition, these VT exit points were in the scar border zone/gray zone as determined by the DE-MRI scans.

Chapter 3. VT Inducibility Prediction: A combined modeling and clinical approach

As local spatial heterogeneity of electrical properties can play an important role in arrhythmogenesis, the differences in AC, APD and ARD-RC between neighbouring regions in the LV (each region ≈ 3mm 3 ) were calculated by evaluating the gradient of the parameters in polar coordinates (as presented on bull's eye plots). The gradient amplitude reflects how smooth or abrupt the transition is between values in the spatial domain. Using a definition of "high gradient" being greater than the 90th percentile of the gradient range, the proportion of LV regions with high gradients was analysed for each patient. There were no statistically significant differences in the number of elements with high gradients (≥ 90 percentile) between patients. However if we look at a combined heterogeneity of both AC and APD-RC in patients, those with inducible VT (P1 and P2) have combined higher values and can then be separated in this representation, see Figure 3.4.

VT Induction and Clinically Observed Exit Points

Patient 1 had induced sustained monomorphic VT with a cycle length (CL) of 275ms at Stage 11 of the Wellens protocol (S1 500ms, with S2, S3 and S4 at 250ms, 240 and 360ms respectively). Patient 2 had induced SMVT with a cycle length of 245ms at Stage 4 of the Wellens protocol (S1 600ms, with S2 at 400 and S3 at 360ms). In Patient 1, the endocardial activation recorded from EAM showed the re-entrant VT circuit initiating from the LV lateral wall, spreading anteriorly and then posteriorly before returning to the lateral wall. In Patient 2, the re-entrant VT circuit exited from the LV infero-lateral wall, spreading antero-laterally and then towards the septum, before returning to the infero-lateral wall. The patterns of activation for both re-entrant VT circuits are illustrated in Figure 3.5.

Analysis of the UEG recording from the EAM showed a higher percentage of maximum APD-RC slope that is 2SD above the mean LV APD-RC slope at the observed exit points (32%) compared to the scar (4%) and non-scar regions (0.2%). The absolute values in APD and tissue conductivity between the three areas were not as distinct although the exit point values (APD: median 319ms, IQR 45ms; AC: median 3.4m 2 /s 2 , IQR 2.5m 2 /s 2 ) appeared to lie between those of the scar (APD: median 385ms, IQR 136ms; AC: median 0.4m 2 /s 2 , IQR 2.1m 2 /s 2 ) and non-scar regions (APD: median 305 IQR 49ms; AC: median 4.7m 2 /s 2 , IQR 2.4m 2 /s 2 ). It is possible that it is the heterogeneity in regional tissue electrical properties distributions rather than the absolute mean or median that may allow us to identify the exit points from other regions of the myocardium. However, this observation is based on just two subjects with induced VT during clinical VT simulation study.

Model-Predicted vs. Clinically-Observed Induced VT

The VT stimulation studies performed on Patients 1-3 were simulated in silico using the personalized computer model. Sustained monomorphic VT was induced at Stage 7 of the virtual Wellens protocol for Patient 1 (drive-train S1 at 400ms, first coupled extra stimulus S2 at 220ms) and stage 5 for Patient 2 (drive-train S1 at 500ms, and cycle length of 260ms compared with the clinically observed cycle length of 275ms with a macroscopically similar activation pattern and a predicted exit point that matched with the observed clinical one. The personalized model of Patient 2 also predicted a positive VT (induced with drive-train S1 at 500ms, S2 at 200ms) stimulation study with a cycle length of 250ms compared with 245ms for the clinical VT. The induced VT was sustained and the re-entrant pathway stable. Notably the direction of the activation pattern during the predicted re-entrant VT was reversed from that observed clinically; however, the predicted exit point correlated with the clinically observed one. We performed additional quantitative analysis on the patients' 3D anatomies to make a comparison of clinically observed and predicted exits points in terms of anatomic location. We defined the exit region compassing anatomical points with an activation time within 10ms of the earliest activation. The distance between the clinically observed and simulated exit points was defined as the 3D Euclidean distance between the centres of both exit regions. A difference of 7mm and 8mm were found, respectively, for patients 1 and 2. The personalized

Simulated VT Stimulation from Additional Sites

VT stimulation study was also performed in silico by pacing from other sites than the RV apex used during the clinical study for Patients 1 and 2. Different VT circuits with three additional exit points were observed in both patients. The exit points were located on the boundary of scar in the region of the gray zones from the DE-MRI images. Similar to the clinically observed exit points, they were mostly in the region of maximum APD-RC slope. The composite of different exit points are plotted in Figure 3.7 in terms of their APD restitution properties and AC.

Three-Dimensional VT Circuit Visualization

The computer model enabled the prediction of a 3D VT circuit, as opposed to the 2D VT activation pattern observed by the non-contact EAM, by taking into account the 3D geometric information from MRI including the transmurality of scar core and gray zone across the LV wall. This allowed additional insights on the wave propagation within the myocardium through computer simulation. For Patient 2, the main VT circuit in the myocardium computed from the geodesic path of the predicted VT isochrones is shown in red in Figure 3.8. Between the entry (the latest activated region) and the exit (the earliest activated region) points, the activation wave front of the re-entrant VT propagates in the region of the gray zone. The estimated wave front path surrounds the scar core and lies within the scar border zone, entering from the endocardial surface, meandering within the ventricular wall, and exiting via the epicardial surface.

Results for solely image-based personalized parameters.

We performed additional sensitivity analysis using solely image-based personalized electrical parameters instead of personalized patient-specific electrical data in the two patients with positive clinical VT stimulation studies and found that using combined personalized electrical and image data can potentially improve the accuracy of VT inducibility and predictions regarding the location of exit zones predictions. Using non-personalised empirical electrical parameters VT was not inducible in Patient 1. It also induced a different macroscopic VT circuit morphology from the one seen during the clinical study in Patient 2. The Euclidean distance between the centre of the clinical exit region and that predicted using non-personalized electrical parameters was 37mm, and with using personalized electrical parameters it was 8mm.

Figure 3.8: Estimation of the intra-myocardial path between the entry and exit points during re-entrant VT. Estimated endocardial 2D geodesic path between the entry and exit points from electroanatomical study during clinical VT stimulation study is demonstrated by the isochrones map. Estimated intra-myocardial 3D geodesic path (red line) takes into account of the scar heterogeneity from the high-resolution MR images. The entry and exit points are illustrated by the graded colours on the 3D path with blue denoting late activation and red denoting early activation. The activation wave path is shown traveling within the region of scar border zone (purple wire mesh), between the scar core and the gray zone.

Discussion

The present study provides new insights into the prediction of VT circuits using a computerized biophysical model. The principal findings of our study are as follows:

1. The location of scar border zone (gray zone) on high resolution DE-MRI correlates with areas of abnormal measured and model derived electrical properties.

2. The model provided mechanistic insights into the underlying electrophysio-logical substrate of the VT circuits showing that the VT exit points have a substantially higher percentage of tissue with steep APD-RC slopes compared to surrounding tissues.

3. In silico VT stimulation studies were able to predict inducibility of VT in patients at risk of ventricular arrhythmia.

4. The characteristics of the modeled VT circuits correlated well with the clinically observed circuits in terms of cycle length, macroscopic activation patterns and VT exit sites.

Co-location of Tissue and EP Properties Heterogeneity

We found a correlation between the spatial heterogeneity in APD, APD-RC slopes and AC and the location of gray zone detected on high-resolution DE-MRI. This suggests that an admixture of scar and myocardium alters the regional restitution properties and provides potential substrate for arrhythmia. In the two ICM patients who developed sustained VT, the VT exit points had 32% of tissue with steep restitution slope amidst areas of tissue with shallow slopes. These two patients also had substantially greater overall spatial heterogeneity in APD and APC-RC slopes. These results suggest that VT exit points are co-located with heterogeneous APD-RC slopes. Ciaccio et al. has demonstrated potential mechanisms using a geometric model that could account for this co-location [Ciaccio 2007]. Furthermore, the presence of clumps of fibroblasts, although electrically non-excitable, results in a slowing of the electrical wave front propagation in the myocardial tissue in which they are embedded. These fibroblasts can modulate cellular ion channel remodelling and therefore tissue electrical properties through a variety of mechanisms including mechano-electrical feedback via stretch-activated ion channels, close coupling of nearby cardiac myocytes via connexin, and by altering the orientation of surrounding myocardial fibres. Some of these mechanisms have been shown in vitro, but remains speculative in vivo studies. However, collectively, these effects may provide a substrate for ventricular arrhythmia [Mitchell 2003] [Fenton 1998] [Camelliti 2005]. Thus, the suggested relationship between the degree of heterogeneity and propensity for sustained ventricular arrhythmia may not be linear, as the absolute differences in electrical parameters expressed as gradients between neighbouring domains in LV did not distinctively identify those patients at risk of sustained ventricular arrhythmia in the present study.

The present study results support recent animal and human clinical studies that have shown potential target sites for successful ablation for VT lie in areas of heterogeneous zones of scar and healthy myocardium identified on DE-MRI [Estner 2011] [Perez-David 2011] [Dickfeld 2011]. These findings are also in keeping with previous studies that showed spatial heterogeneity in APD and APD-RC slope within the LV [Yue 2005[START_REF] Nash | Whole heart action potential duration restitution properties in cardiac patients: a combined clinical and modelling study[END_REF]]. These studies did not have the benefit of high resolution MRI scar information, and so could not explain the anatomical cause of this heterogeneity in electrical properties. Nash et al. speculated that regional stretch strain pattern resulting from ischemia might have led to electrical remodelling through mechano-electrical feedback. As in the two patients they presented in detail, the regions with the greatest spatial heterogeneity corresponded to the territories of the diseased coronary arteries. It may well have been that these regions had areas of gray zone. The mean APD and APD-RC slope values were greater in our cohort, compared to the other studies which showed a skew towards lower APD-RC slope values [Yue 2005[START_REF] Nash | Whole heart action potential duration restitution properties in cardiac patients: a combined clinical and modelling study[END_REF]]. One explanation is that our study cohort had advanced cardiomyopathy with adverse remodelling, illustrated by the severely impaired LVEF. Computation modelling and experimental studies have demonstrated that theoretically APD-RC slope >1 is needed to result in electrical instability and initiate ventricular fibrillation [Nolasco 1968[START_REF] Karma | Electrical alternans and spiral wave breakup in cardiac tissue[END_REF][START_REF] Gilmour | [END_REF]] [Garfinkel 2000]. Others have demonstrated that it is the heterogeneity in APD-RC slopes that is needed for initiation and stability of the re-entry mechanism [START_REF] Nash | Whole heart action potential duration restitution properties in cardiac patients: a combined clinical and modelling study[END_REF]] [Clayton 2002]. These data lend support to our observed functional heterogeneity in myocardial tissue restitution properties at the VT exit points.

In silico VT stimulation studies in patients

Encouragingly, the personalized cardiac models encompassing both anatomical and electrical properties in the present study were able to predict not only the inducibility of re-entrant VT, but also the re-entrant VT circuit properties and anatomical locations of the substrate. Whilst cardiac modelling has been an active research area for decades, personalized cardiac modelling using patient-specific clinical data is in its infancy. Recently, Ashikaga et al. presented a retrospective study that compared the actual target VT ablation sites with that predicted from an image based cardiac model and found a good correlation between the predicted ablation sites located mostly at the scar border zone and the actual ablation sites [Ashikaga 2013]. In keeping with our data, they highlighted an advantage of incorporating 3D geometric information gained from MRI, which allowed the prediction of the potential critical isthmus on the epicardial surface. This is important, as conventional mapping techniques used during ablative therapy are limited especially when the substrate of the VT circuit lies on the opposite side of the mapping surface. Our study is unique in that as well as anatomical and scar information, we incorporated functional data including APD-RC properties and tissue conductivity and their spatial heterogeneity into our patient-specific models. Previously, Arevalo H et al. suggested that VT dynamics were primarily governed by the geometric parameters of the scarcore and border zone using image-based computational VT modeling work with scar data acquired from high-resolution MRI of ex vivo canine heart (reconstructed voxel size 200µm 3 ) and empirical ionic cell parameters [Arevalo 2013]. Whilst we realize that contrast-enhanced MRI provides important scar geometry that governs the substrate of re-entrant VT, we also recognize that the resolution of the current standard clinical MR imaging technique is limited (common voxel size 2×2×8 mm 3 ) clinical approach in providing the level of spatial geometry details that we would like to see at the border zone. Thus the term "gray zone" reflecting the spectrum of intermediate signal intensity seen on MRI scar images at the scar border zone was originated. We believe that additional knowledge and understanding of patient-specific heterogeneity in local electrical parameters would assist in predicting the likely culprit conduction channels/isthmus, and not the bystanders, that is critical to the clinical VT. This by confirmed by the results of our simulations on the solely image-based personalization of the electrical parameters in the model. We found that using combined personalized electrical and image data can potentially improve the accuracy of VT inducibility and predictions regarding the location of exit zones predictions.

Clinical application

Potential for Circuit Prediction using Personalized Computer Models to Guide Ablation. Successful VT termination through ablation may be achieved when the critical isthmus is successfully interrupted with ablation lesions. Predicted VT exit points could be potential targets for ablation. In silico personalized models may offer significant clinical benefit in predicting the risk of ventricular arrhythmia in patients and guiding treatments including ICD implantation and VT ablation. If VT exit points can be predicted with biophysical models then this information may be used to guide ablation in patients. Such models also offer additional flexibility as the model can simulate any combination of paced stimuli from different locations with varying pacing cycle lengths, which may not be feasible in clinical practice. We have indeed demonstrated that additional potential exit points, which still lay at the scar border zone/gray zone, can be induced from in silico VT stimulation studies.

Study Limitations

The study is limited by the small number of patients included, with only a subset of patients undergoing a VT stimulation procedure. However the invasive nature of the study precludes analysis of a large number of patients. The electro-anatomical data are derived from non-contact mapping with the inherent limitations of this type of mapping [Schilling 1998]. Approximately 13% of the EGM recording from each non-contact mapping with indiscernible T waves were not suitable for analysis. We used the "alternative method" as oppose to the "Wyatt method" for estimation of APD because we found it was more suited for the virtual unipolar EGM derived from non-contact array mapping as others have found [START_REF] Arthur M Yue | Determination of Human Ventricular Repolarization by Noncontact Mapping Validation With Monophasic Action Potential Recordings[END_REF]]. Though there are additional data in support of validity of the Wyatt method [Wyatt 1981]. Non-contact mapping was chosen for its ability to provide 'beat-to-beat' mapping in the setting of rapid and unstable VT circuits, which could not be mapped with contact mapping techniques. Non-invasive body surface mapping may be incorporated to routine simple electrophysiology study to gain such personalized whole heart electroanatomic data in order to facilitate the translation of the biophysical cardiac model processing pipeline to clinical practice [START_REF] Rudy | Noninvasive electrocardiographic imaging of arrhythmogenic substrates in humans[END_REF][START_REF] Zettinig | [END_REF]]. In addition, estimation of APD-RC properties from single pacing and from multiple recording sites with varying distances has a degree of limitation as the measurement from the distal site has some attenuations due to the premature stimuli at varying short coupling intervals impacting on conduction velocity restitution. However, this confers smooth systematic influences across the LV myocardium which is similar in all patients studied. The model required high-resolution MRI data. In patients with pre-existing ICDs such imaging may not always be feasible. The cardiac model also made several simplifications in particular for the Purkinje network and pathological changes of cardiac fibre orientations. Though lacking detailed representation of such microscopic parameters, the model exhibited the main macroscopic properties of the tissue such as conductivity, APD and APD restitution. As with any personalization of computational physiology model, there is a degree of uncertainty due to the limitations in acquired temporal and spatial resolution clinical data, as well as in the model. Computational methods were developed to incorporate a simpler model in an efficient Bayesian inference method that can take into account the uncertainties on the data and model parameters in the personalization [START_REF] Konukoglu | [END_REF]]. Extending such methods to the model used in the study could account the uncertainties in the application of the model.

Conclusion

Patient-specific spatial heterogeneity of restitution properties were the distinguishing features of ventricular arrhythmogeneicity, with re-entrant VT exit points present in regions of higher maximum APD-RC slopes compared with the surrounding tissue. These regions were within the gray zone identified by DE-MRI. Our personalized biophysical model was able to predict macroscopic VT circuits and exit point locations in agreement with clinically observed datasets. Our results suggest that patient-specific cardiac models may offer incremental clinical benefit in terms of ventricular arrhythmia risk stratification and in the planning and delivery of ablation strategies for re-entrant ventricular arrhythmias.

Introduction

Sudden cardiac arrest (SCA) is a leading cause of death in the world, with 350,000 deaths per year in the USA, and similarly in Europe. Its main cause is cardiac arrhythmia with RFA increasingly being used to treat it but with an unsatisfying success rate due to the difficulty to find the ablation targets. There is therefore a need to identify the arrhythmia substrates and the optimal ablation strategy to substantially improve its success rate.

Animal myocardial infarction models have been used to explore the poor coupling of surviving cells within fibrotic scar to the rest of the myocardium [START_REF] Reddy | [END_REF]]. They are known to generate regions with critically slow conduction properties and to promote electrical circuit re-entry. Similarly, in humans, inhomogeneous scars and/or heterogeneous necrosis seen in ischemic and dilated cardiomyopathy provide an ideal substrate for ventricular tachycardia (VT) [START_REF] Jaïs | [END_REF]]. These bundles of surviving tissue can be identified using electrophysiological (EP) mapping, a lengthy and invasive method that records cardiac electrical signals through intra-cardiac catheters.

Currently, clinicians rely in the identification of a re-entry circuit through activation and entrainment mapping to ablate critical isthmus. Nonetheless, a considerable amount of patients is non-inducible or have poor hemodynamic tolerance for VT mapping. The clinical study in [START_REF] Jaïs | [END_REF]] describes a distinct approach where the isolation or dissociation of electrically surviving myocardial fibers within scar regions is performed through the ablation of local abnormal ventricular activities (LAVA) recorded during sinus rhythm. LAVA were defined as sharp fractionated bipolar potentials occurring during or after the far-field electrogram and which arise from pathological tissue. Its complete elimination through RFA has been associated with a superior clinical outcome [START_REF] Jaïs | [END_REF]].

Delayed-enhanced magnetic resonance imaging (DE-MRI) enables a non-invasive 3D assessment of scar topology and heterogeneity with millimetric spatial resolution. It has been hypothesised that areas of intermediate signal intensity in LE-MRI, referred to as the grey zone, are likely to host both scarred and surviving myocardium related to arrhythmia in ischemic populations [START_REF] Jaïs | [END_REF]] [Perez-David 2011]. Some studies [Fernández-Armenta 2013] have focused on the detection and characterization (length, orientation, etc.) of grey zone channels with surviving fibers that are surrounded by scar core with the use of high-resolution DE-MRI. The study in [START_REF] Sasaki | [END_REF]] sought to quantify the association of scar on LE-MRI with local electrograms and vetricular tachycardia circuit sites in patients with non-ischemic cardiomyopathy. Image features such as scar transmurality and intramural scar type were associated with electrogram characteristics, including amplitude and number of deflections. However, these feature sets fail to account for complex texture patterns that can can arise due to the intertwining of healthy and scarred tissue than can be responsible for these electrophysiological pheonomena.

Machine learning techniques have successfully found applications in medical image processing and have since played an essential role in the field [Criminisi 2011a] 4.2. Clinical Data 53 [START_REF] Lempitsky | Random forest classification for automatic delineation of myocardium in real-time 3D echocardiography[END_REF]]. In this work, they will aid in relating imaging features to LAVA occurrences with direct clinical applications to pre-emptively define mapping and ablation targets, to increase the success rate and to decrease the duration of such procedures (currently>6h).

Specifically, we evaluated the predictive power of locally computed intensity and texture-based MRI features to identify RFA targets. On the methodological side, we used random forests with advanced image features and classifier parameter estimation using nested cross-validation. However, using machine learning approaches on such complex multimodal data can prove difficult due to the inherent errors in the training set. We present a detailed analysis of these error sources, due in particular to catheter motion and the data fusion process, their formal integration in the training process, which is rarely done in machine learning approaches, and demonstrate an improved algorithm performance.

Clinical Data

Three patients referred for cardiac ablation for post-infarction in ventricular tachycardia were included in this study. They underwent cardiac MRI prior to highdensity EP contact mapping of the endocardium (Patients 1 and 2) and epicardium (Patient 3). Signals were categorised as normal or LAVA by an experienced electrophysiologist following the definition of LAVA signals detailed in [START_REF] Jaïs | [END_REF]] and [Sacher 2014]. Figure 4.2 shows characteristic electrical readings from two electroanatomical points (EAP) depicting healthy and a LAVA electrograms. Our inhome software facilitated the visualization of all electrode readings from the CARTO study and allowed the electrophysiologist to review the signal labeling and correct for any inaccuracies. In these images, the top row depicts the bipolar signal obtained from the difference of the two distal-most electrodes in the catheter (M1-M2) and it is the one commonly used for the identification of LAVAs. Nonetheless, unipolar signals and other bipolar measurements can be chosen to be displayed complement electrical behaviour understanding and to clarify difficult cases for labeling. The volume of tissue that influences the recording at a particular EAP is referred to as the catheter's sensing range and in this study it is represented by a sphere of 10mm radius, illustrated in Figure 4.3. Because of this phenomenon, two sufficiently close electrical readings would be influenced by roughly the same tissue section in the myocardium. Additionally, because of uncertainties which will be detailed in Section III, a spatial filtering of the EP study was made. Under these hypotheses, if two spatially contiguous (less than the catheter's sensing range) EAPs were found to have different labels (LAVA and healthy), the EAP assigned the healthy label would be removed from the study. This would be justified in a clinical scenario as the objective in a LAVA driven ablation procedure is to identify and ablate the totality of regions presenting LAVA characteristics. Table 4.1 summarizes the characteristics of the electrophysiological datasets for each of the patients in this study.

Electrophysiological (EP) Data

Imaging Data

Image Acquisition: The scar was imaged on a 1.5 Tesla clinical device (Avanto, Siemens Medical Systems) 15 minutes after the injection of a gadolinium contrast agent. A whole heart image was acquired using an inversion-recovery prepared, ECG-gated, respiratory-navigated, 3D gradient-echo pulse sequence with fatsaturation (1.25×1.25×2.5mm 3 ).

Segmentation:

The myocardium was manually segmented on reformatted images of isotropic voxel size (0.625mm 3 ). Abnormal myocardium (dense scar and grey zone areas) was segmented using adaptive thresholding of the histogram, with a cut-off at 35% of maximal signal intensity. Segmentations were reviewed by an experienced radiologist, with the option of manual correction.

Sources of Uncertainty

We aim to identify differences in regional image characteristics between LAVAinducing and healthy tissue. Nevertheless, despite the integrated catheter localisation and previous registration with anatomical data, two main sources of uncertainty between the electrophysiological measurements and the imaging data remain: data fusion and temporal displacement.

Data Fusion

Meshes generated by EP mapping systems are a rough approximation of the shape of the ventricular cavity, as seen in Figures 4.4 and 4.5. Therefore, a registration is needed between the EP recording locations and the image segmentation. Some works, like [Reddy 2004], use iterative closest point to register the EP recording locations to the endocardial surface and surface to the aorta to resolve amiguities in rotation of the left ventricle. In our case, the aortic readings were not available. Therefore, this is done manually by using landmarks and matching between low voltage areas and scars. Then the EAPs are projected on the image-based mesh by finding the closest cell, and an evaluation of the registration uncertainty is present in the resulting projection distance (PD), which is given by the Euclidean distance between the EAP and its closest cell center. An illustration of the projection distance maps is shown in on the top panel of

Temporal Displacement

Due to breathing and cardiac motions, the recording catheter is displaced throughout the 2.5 seconds of recording time. Magnitudes varying significantly among electro-anatomical points (EAPs) as is shown in Figure 4.5, where the red ellipsoids describe the catheter's temporal position variation. A map of the distribution of the magnitudes of these displacements projected on the endocardial surface can be found in the bottom panel of Figure 4.6.

Image Feature Computation

The previously described properties allow us to define the scale at which the image features will be looked at and how to quantify the errors.

Intensity-based Features

The work in [Perez- David 2011] showed that MR image intensities can be used to discriminate the heterogeneous substrate. Voxels contained inside the sensing range were used to compute intensity-based features, including minimal, maximal, mean and standard deviation values. Another feature, defined as the standard deviation over the average intensity in the region, was included. Myocardium thickness was calculated and the scar transmurality was defined as the extent of scar through the entire myocardium thickness.

Texture-based Features

Texture from medical images provide us with information inherent to the underlying tissue. This information can later be used in applications as cardiac image indexing and retrieval [START_REF] Glatard | Texture based medical image indexing and retrieval: application to cardiac imaging[END_REF]], or to predict vascular events [van Engelen 2014]. Grey level co-occurence matrices (GLCM) are matrices of the joint probability of occurrence of a pair of grey values separated by a displacement d = (dx, dy, dz). Haralick features are statistics computed on GLCM that emphasize specific texture properties and have been extensively used in medical image analysis [Tesar 2007]. In our study, the GLCM were computed around the center of the myocardium were the EAP had been projected using a ROI of window size of 11×11×11 pixels (∼9.4×9.4×9.4mm). Three distances from the central pixel (1, 2 and 4 pixels), 13 directions and 12 Haralick features were considered, resulting in a 468 element texture feature vector per EAP analyzed.

Concatenation of the intensity and texture features yielded a final image-based feature vector of 475 dimensions which was used for classification.

Confidence-based Learning Framework

Random Forests Classifier

Random forests are discriminative classifiers created in an intuitive and easily understandable structure that also provide informative uncertainty measures on the classification results [Criminisi 2011a]. They have successfully found multiple applications in medical image processing [Criminisi 2011a] [ [START_REF] Lempitsky | Random forest classification for automatic delineation of myocardium in real-time 3D echocardiography[END_REF]] [Zikic 2013]. We used the Python implementation from the scikit-learn library [Pedregosa 2011].

Nested Cross-validation

Cross-validation has been shown to be among the best ways to estimate performance [Kohavi 1995]. It predicts the fit of a model to a hypothetical validation set when an explicit one is not available. Also, because all observations are used for both training and validation (though never within the same fold), and each observation is used for validation exactly once, this method exploits to the fullest the information in a set with a low risk of overfitting, particularly important in small datasets like ours. We used stratified cross-validation and optimized the classifier for precision performance. The use of nested cross-validation, with a parameter-tuning inner loop and an outer loop for performance estimation, avoided an optimistic bias introduction into generalization estimate [Ruschhaupt 2004].

Feature Selection

Due to the high dimensionality of our feature vector, the effect of feature space reduction was next assessed. Univariate t-Test statistics were used to assess feature significance [Guyon 2003]. Three reduced datasets were created, including the 50%, 4.5. Confidence-based Learning Framework 59 25% and the 10 most relevant features (MRF). Additionally, a feature subset was generated containing only the intensity-based features.

Uncertainty Assessment

Uncertainty assessment has become a crucial aspect in our field, particularly to account for the impact of noisy or inaccurate ground truth when being fed into learning algorithms. It arises in the form of weights assigned to the contribution of multiple sources to perform anatomical segmentations [van Opbroek 2015] to the impact of the uncertainty in these segmentations on the estimation of physiological parameters [Sankaran 2014].

We derived a principled analysis of confidence impact on classification. Inspired by cost-sensitive learning, we formulate the problem as samples (x, y, c) drawn from a distribution D on a domain X × Y × C with X being the input feature space, Y corresponding to the binary output class and C to the confidence associated with each sample. We aim to learn a classifier h : X → Y which minimizes the new expected classification error:

E(x, y, c ∼ D)[cI(h(x) = y)] (4.1)
Using the Translation Theorem 2.1 in [Zadrozny 2003] we can compute and draw samples from a distribution D ′ such that the optimal error rate classifiers for D ′ are optimal cost minimizers for data drawn from D. We derive how this modifies the training using weights to simulate the expectation of finite data E(x, y ∼ D)[f (x, y)] as:

E(x, y ∼ D)[f (x, y)] = 1 c cf (x, y) (4.2)
equivalent to importance sampling for D ′ using distribution D, so the modified expectation is an unbiased Monte Carlo estimate of the expectation with respect to D ′ [Zadrozny 2003]. In random forests, the node split criterion is information gain:

IG = H(S) - i=1,2 |S i | |S| H(S i ) (4.3)
with |S| being the number of samples in a node before split, |S i | being the number of samples of each children node and H(S) the Shannon entropy:

H(S) = - c∈C p(c)log(p(c)) (4.4)
where p(c) is calculated as normalized empirical histogram of labels corresponding to the training points in S, p(c)

= |S i | |S| .
Using weighted instances, p(c) is replaced by p w (c), which has the following formulation:

p i w (c) =

Weights of samples of class c in node i

Weights of samples in node i (4.5)

p i w (c) = Sc W i S W i (4.6)
This yields a sample weighted formulation of the information gain that can be written as:

IG = H(W ) - i=1,2 S W i S W H(W i ) (4.7)
where W are sample weights at the parent node and W i are sample weights that have been passed to each child node. H(W ) is given by:

H(W ) = - c∈C p w (c)log(p w (c)) (4.8)
To our knowledge, it is the first time such formulation is derived in this context and we believe it strengthens our approach's methodological ground.

We analysed the influence of two factors affecting the certainty in EAPs and imaging data correspondences by weighting our training samples according to the confidence we have on their image features.

Projection Distance. More confidence is assigned to the imaging features computed from EAPs with low PD with respect to those with high values. It is defined as the Euclidean distance between the EAP and the center of the cell in the mesh closest to the given point:

P D = CellCenter -EAP (4.9)
Temporal Displacement. The covariance of the position matrix is obtained and an ellipsoid with radii 2 diag(D) is generated, where D is the matrix containing the eigenvalues along the main diagonal. The major ellipsoid radius defines the temporal displacement, as following:

T D = max(2 diag(D)) (4.10)
Intuitively, image features from EAPs with smaller major ellipsoid radius are more reliable as they are less affected by movement.

Each EAP is assigned a confidence value by linearly scaling either the PD or the temporal displacement to a weight parameter with range of [0.5, 1] where 0.5 corresponds to the lowest confidence and 1 to the highest. Additionally, a combined uncertainty weight is defined as the product of both uncertainty sources. We explore the uncertainty inherently introduced to our dataset due to these factors with three extra experiments by using the previously described sample weighting schemes.

Evaluation Metrics

The results for the classification results will be assessed using precision-recall (PR) and receiver operating characteristic (ROC) curves. A PR curve illustrates the trade-off between the proportion of positively labelled examples that are truly positive (precision) as a function of the proportion of correctly classified positives (recall). A good performance line lies in the upper-right portion of the graph and has a high area under the curve (AUC) value.

A ROC curve depicts relative trade-offs between benefits (true positives) and costs (false positives). Any deviation from the diagonal (representing random guess and with 0.5 AUC) into the plot's upper triangle represents a better performance than chance. Classifiers with higher ROC AUC values are said to have a better average performance.

Results and Discussion

The results of classification using the full feature set for each patient are shown in By closely looking at the PR curves, we can conclude that the classifier is able to retrieve approximately more than half of the LAVA instances without having a considerable drop in the precision. It is when the totality of the LAVA instances are recovered that precision is compromised. This might correspond to areas in which the LAVA regions are spatially close to healthy ones, therefore finer descriptors of the adjacent areas should be explored. Nonetheless, for our approach some precision can be compromised as the ablation procedure might not be able to distinguish between closely spaced cardiac regions. Plots for Patient 3 show that it presents great variability between folds. Some have perfect scores, probably representing the typical LAVA image signatures, while others perform poorly, possibly due to outlier LAVA present in the testing phase. This can also be due to the patient having an epicardial study compared to the endocardial mapping of Patients 1 and 2.

Results Using Subsets of MRF

The classification results obtained using only 50%, 25% and the top 10 MRF in each dataset are summarized in Table 4.2. The purpose of this task was to investigate the redundancy in the feature set. In general, the use of 50% of MRF resulted only in a small drop in AUC PR and AUC ROC scores, suggesting some feature redundancy. Nonetheless, while Patient 2 and 3 had a slight increase in performance when using its 25% MRF, all patients suffered a considerable drop in score when using 10 MRF. Due to the characteristics of our current feature selection scheme (univariate filtering method), it failed to assess groups of features that work together to better discriminate between classes.

The results for the assessment of intensity feature importance w.r.t. the full set of features are included in Table 4.2. Using only intensity features for classification led to a large drop in overall classification performance. This shows that advanced texture patterns are required to describe the complex intertwining of myocardial fibres in scarred and grey zone areas responsible of LAVA generation.

Uncertainty Impact on the Prediction

Results of classification with the full feature set and weighted samples are shown in Table 4.3 and Figure 4.9. A general increase in performance is observed when weighting samples according to their proximity to the location of image feature computation and their temporal position stability. This confirms our hypothesis that a lower confidence should be assigned during training to EAPs with large PD or temporal displacements. Weighting samples with a combination of both uncertainties results in a similar improvement in the classification performance as when the elements were used independently. Currently, this combined uncertainty is a naive product of both elements. A different fusion should be explored to better exploit both uncertainty sources and construct a more reliable estimate of the confidence of a given EAP.

The optimal random forests construction parameters found during nested cross validation with precision optimization and and temporal displacement weighting are included in Table 4 4.5, when depicted as false negative regions as in Figure 4.10 the numbers dropped to 4, 2 and 3, for each of the patients.

Conclusion

We presented the use of a machine learning framework, using intensity and texturebased local imaging features in the vicinity of myocardial scar and grey zones towards the prediction of RFA target localisation. Additionally, we detailed the uncertainty in the data and explored its impact on the classification results. For both PR AUC and ROC AUC, we scored above 0.75 and an extra 0.05 was gained when using uncertainty evaluation to weight the training data. Finally, a preliminary output with visual interpretation and potential use in a clinical environment was presented.

The aim of the current work was to analyse the feasibility of classification using solely image-based features on complex multi-modal data. So far, the choice of features produced encouraging results. We are aware that an increase in the size of the patient database is required in order to aim for inter-patient analysis, but this work serves as a proof of concept with results good and encouraging enough to warrant further investigation and open up possibilities for non-invasive cardiac arrhythmia ablation planning.

Extension to other pathologies

A similar learning algorithm was applied in the collaboration mentioned in Appendix A to differentiate between rabbits presenting LQT2, E4031 or LMC conditions based on image-based mechanical data.

Introduction

Invasive measures of cardiac activity obtained during an electrophysiology (EP) study provide insightful information about the electrical characteristics of the analysed myocardium. Under pathological conditions, these extracellular recordings deviate from the basic (healthy) signal shape and present multiple deflections and fractionation episodes.

More than a decade ago, fractionated electrograms (EGMs) were thought to be caused mainly by artifacts related to the electronics of the acquisiton system. Although artifacts indeed may cause complex electrograms, most of the complex and fractionated electrograms are caused by the peculiar behaviour of activation fronts, due to structural and electric complexity of the underlying tissue [de Bakker 2010]. Previous studies have used synthetically generated EGMs to explore different pathological phenomena. The authors in [Ciaccio 2014] describe the generation of electrogram fractionation from changes in activation wavefront curvature in experimental canine infarction. In [Álvarez 2012], the modelling of intracardiac recordings was used to aid in the reconstruction of cardiac ischemia. [START_REF] Keller | [END_REF]] studied the influence of different catheter angles, locations and filter settings on the morphology of simulated intracardiac electrograms and compared them to clinical signals.

On structurally diseased hearts with fibrotic scar, bundles of surviving tissue promote electrical circuit re-entry and are a cause of arrhythmias. Local abnormal ventricular activities (LAVA), sharp fractionated bipolar potentials occurring during or after the far-field electrogram, have been shown to indicate surviving fibres within the scar and have been successfully used as targets for radio-frequency ablation [START_REF] Jaïs | [END_REF]].

In this work, we test the feasibility of using DE-MRI image-based simulation to reproduce LAVA-like patterns in catheter recordings of intracardiac electrograms. Through this work, we aim to shine some light on the understanding of the mechanisms underlying electrogram fractionation in the grey zone.

Clinical Data

Five patients referred for cardiac ablation for post-infarction in ventricular tachycardia were included in this study. They underwent cardiac MRI prior to high-density EP contact mapping of the endocardium.

Electrophysiological Data

The CARTO mapping system (Biosense Webster) enables the 3D localization of the catheter tip and provides the distribution of EP signals on cardiac surfaces. Contact mapping was achieved in sinus rhythm on the endocardium (trans-septal approach) with a multi-spline catheter (PentaRay, Biosense Webster). Signals were categorised as normal or LAVA by an experienced electrophysiologist.

Imaging Data

The scar tissue was imaged on a 1.5 Tesla clinical device (Avanto, Siemens Medical Systems) 15 minutes after the injection of a gadolinium contrast agent. A whole heart image was acquired using an inversion-recovery prepared, ECGgated, respiratory-navigated, 3D gradient-echo pulse sequence with fat-saturation (1.25×1.25×2.5mm 3 ).

The myocardium was manually segmented on reformatted images of isotropic voxel size (0.625mm 3 ). Abnormal myocardium (dense scar and grey zone areas) was segmented using adaptive thresholding of the histogram, with a cut-off at 35% of maximal signal intensity. Segmentations were reviewed by an experienced radiologist, with the option of manual correction.

Cardiac Electrophysiology (EP) Simulation

As described in Section 2.5, the cardiac electrophysiology modeling field proposes several mathematical tools to represent the electrical activity in the heart. For the purposes of this study, we chose a model that was able to represent complex cardiac electrical phenomena while keeping the number and variation of the involved variables tractable. We therefore chose the Mitchell-Schaeffer [Mitchell 2003] model for EP simulation. Its main advantages include: (1) its simplicity, as it only includes two differential equations; (2) the relation between its parameters and physiological behaviour make it easier to understand and interpret; (3) its ability to simulate arrhythmia macroscopically due to its restitution parameters.

The chosen model has two variables: u, the transmembrane potential and z, a secondary variable in charge of controlling the repolarisation phase. The model is governed by the following two equations:

       ∂ t u = div(D∇u) + zu 2 (1-u) τ in -u τout + J stim (t) ∂ t z = (1-z) τopen if u < u gate -z τ close if u > u gate (5.1) Where D = d • diag(1, r, r
) is an anisotropic diffusion tensor that enables the conduction velocity in the fiber direction to be 2.5 times greater than in the transverse plane R = 1 2.5 2 . The parameters τ in and τ out define the repolarisation phase whereas τ open and τ close are responsible for the gate opening or closing depending on the change-over voltage u gate . J stim is the stimulation current at the pacing location.

The SOFA public framework has a GPU implementation of the Mitchell-Schaeffer EP model and it was chosen as the simulation framework for this Chapter 5. Image-based Simulation of LAVA Intracardiac Electrograms work. Technical details on the implementation are more thoroughly described in [START_REF] Talbot | [END_REF]].

Pseudo-personalisation: Zone-specific EP Properties

The electrophysiological properties of the infarcted and border zone tissue are known to differ from those of the healthy myocardium. For the construction of our personalized model, DE-MRI was used to asses the viability of the cardiac tissue and to identify tissue as either healthy myocardium, scar or grey zone. Segmentations were done by an experienced radiologist. The derived volumetric ventricular model incorporating this information is shown in Figure 5.1. Conductivity in the grey zone was decreased by 90% with respect to its value in the healthy myocardium (0.1 × d healthy ) [Arevalo 2013]. Studies in infarcted hearts [Pu 1997[START_REF] Dun | [END_REF][START_REF] Jiang | [END_REF]] [Decker 2010] have reported variations in ionic currents in the border zone tissue generating action potentials that differ from the healthy myocardium in the following aspects:

• 32% lower peak action potential amplitudes • 31% smaller upstroke velocity • 25% longer action potential duration (APD)

In order to achieve the following effects in our model, the term J in = zu 2 (1-u) τ in in Equation 5.1 was modified to zu 2 (a-u) τ in , where variable a controls the peak amplitude for the action potential. Also, to obtain a smaller upstroke velocity, the time constant in the grey zone was set to be τ in GZ = 1.31 × τ in healthy .

The analytical expression for the APD for the Mitchell-Schaeffer model is given by:

AP D = τ close ln τ out 4τ in (5.2)
Implying a linear relationship between parameter τ close and the APD, therefore the values for this parameter in the grey zone were increased by 25%. Finally, because τ in had been modified to account for a smaller upstroke velocity, τ out was modified in the same proportion to keep the ratio τ out /τ in constant.

Scar tissue was modeled as having a null reaction component on the Mitchell-Schaeffer model by setting much higher values of parameters τ in and τ out with respect those of the healthy myocardium. Nonetheless, the diffusion component was left with the values of the grey zone as has been previously described.

Healthy myocardial values were set to the default values appearing in [Mitchell 2003]. Table 5.1 shows the resulting set of parameters used for our personalized simulations.

Sinus Rhythm Computation from Surface ECGs

The QRS window was detected using all six surface electrocardiograms (ECGs): limb leads I, II and III and the augmented limb leads aVL, aVR and aVF, the process is shown in Figure 5.2. A third order bandpass digital Butterworth filter with a bandpass frequency range of 0.5 to 45 Hz was applied to all signals to filter out noise. This filter was chosen because it has a flat frequency response in the passband, the reason for which it is also referred to as a maximally flat magnitude filter.

Next, the gradient for each of the surface ECGs was computed. For each individual lead, a candidate QRS window segment was detected where the absolute gradient value was above the 90th percentile of the absolute gradient throughout the entire signal. Because maximum and minimum peak values in the electrocardiogram would correspond to zero crossings in the gradient signal, and that these could pertain to QRS peaks that should be included in the detection of candidate QRS windows, a gap filling process followed. In healthy subjects, a normal QRS segment varies from 70 to 100ms in duration. To group successive peaks within a single QRS complex together, a maximum value of gap width of 10ms was allowed. Finally, candidate QRS window segments with a width of less than 20ms were filtered out.

Because a QRS segment appears at the same time in all surface ECGs, a majority voting scheme was introduced to detect the final QRS window segments and to filter out undesirable detections due noise or artifacts in individual leads. All six surface leads had an equal weight and only time points where at least three surface leads agreed in the presence of a QRS window segment were finally labeled as such. Figure 5.2 shows the process of QRS window detection from surface ECGs. The sinus rhythm was extracted by analyzing the number of QRS window segments detected during the 2.5 seconds of recording time and computing the average period between them. For the patient shown in 5.2, the mean heart period was calculated to be 0.83 seconds, corresponding to a sinus rhythm of 72 bpm.

Intracardiac EGM Simulation

Extracellular electrograms arise due to trans-membrane currents occurring from differences in the axial voltage gradient at the interface between activated and inactivated myocardial cells [de Bakker 2010], as is shown in Figure 5.3. In this example, an activation wave front is traveling from left to right throughout the myocardium. At the boundary of the propagation front, the depolarised cell has an intracellular potential of +20mV whereas that of the cell in resting state remains at -90mV, therefore creating a current flow from the activated to the inactivated myocyte. To comply with the preservation of charge, the back of the activation front will present a current flow from the extracellular space into the intracellular space, whereas the current flow direction will be reversed in the regions ahead of the activation front, i.e. a current flow from the intracellular space, through the membrane and into the extracellular space. This flow of current in the extracellular space will generate an extracellular voltage difference and therefore a potential dipole.

Figure 5.3: Schematic showing the generated dipole in the extracellular space as a result of the potential gradients in the intracellular space located at the depolarisation wave front. Consequently, the potential field is positive ahead of the activation front and negative at the back.

The Dipole Approach

The notion of resulting dipoles at the depolarisation wave front can be used in order to compute the resulting potential at a given coordinate in the extracellular space. This approach had already been used by [Berenfeld 1998] to compute pseudo pre-cordial and limb-lead ECGs from a Purkinje muscle model. In this study, we introduce a similar approach to simulate unipolar and bipolar intracardiac EGMs at a given position representing the catheter location.

The monodomain formulation considers that the extracellular potential is grounded and thus the intracellular potential is equal to the membrane potential [Clayton 2004]. It is governed by:

C m ∂v ∂t + I ion = ∇ • σ∇v (5.3)
where v represents the transmembrane potential, C m is the membrane capacitance, σ corresponds to the local conductivity and I ion is the current through the cell membrane per unit of area. We define the equivalent current density j eq as: Chapter 5. Image-based Simulation of LAVA Intracardiac Electrograms j eq = -σ i ∇v

(5.4)

-j eq behaves like a flow source density, and it can also be seen as a dipole moment per unit of volume: p = j eq dV (5.5)

According to the volume conductor theory [START_REF] Malmivuo | Bioelectromagnetism: principles and applications of bioelectric and biomagnetic fields[END_REF]], the electric potential registered at a distance R in a homogeneous volume conductor of conductivity σ electrode outside the region occupied by the volume source V is :

Ψ(R) = 1 4πσ electrode V j eq • ∇( 1 R )dV (5.6)
In our case, we model the moving propagation front as a local dipole (Figure 5.3). We first assume that the dipole momentum p X at a position X is spatially invariant over a single tetrahedron H :

∀X ∈ H, p X = p H (5.7)
Second, we assume that we can take into account all the myocardial regions, the ones being activated, the ones depolarised and the ones at rest. The regions that are not activated give almost null dipole momentums and it simplifies the overall calculation.

We discretize Equation (5.5) in space (H tetrahedra of the myocardial mesh) and get the following formulation of the dipole moment of the charge in the volume V H of tetrahedron H :

p H = V H j t eq,H = V H σ H ∇v H (5.8)
The gradient of the electric potential ∇v(X H,center ) at the tetrahedron center X H,center is estimated using the node positions X H,i and the shape vectors D i of the tetrahedron H defined as:

D i = s V H ((X H,i⊕2 -X H,i⊕1 ) ∧ (X H,i⊕3 -X H,i⊕1 ) (5.9)
where s = 1 for i = 2, 4, else s = -1 (Figure 5.4). These shape vectors are useful for evaluating spatial derivatives [Sermesant 2003]. The gradient of the electric potential writes as:

∇v H = ∇v(X H,center ) = 4 i=1 v(X H,i )D i
(5.10)

From Equation 5.6, the contribution Ψ H (X electrode ) of the tetrahedron H to the potential field calculated at a point X electrode is estimated by : Figure 5.4: Schematic showing tetrahedron H, its center X H,center and its nodes X H,i .

Ψ H (X electrode ) = 1 4πσ electrode (V H σ i,H ∇v H ) • (X electrode -X H,center )
X electrode -X H,center 3

(5.11)

Finally, we sum over the whole mesh to get the potential field at X electrode :

Ψ(X electrode ) = N b tetra H=1 Ψ H (X electrode )
(5.12) Bipolar electrograms, on the other side, correspond to the difference in potential between two unipolar measurements. This type of electrograms are useful to study the local activities. The far-field signal is assumed to be similar for both unipolar recordings, therefore, it is largely filtered out [Stevenson 2005].

Unipolar and Bipolar EGM Computations

The recording catheter used in the clinical environment has four electrodes named M1, M2, M3 and M4 from the distal to the proximal. Two bipolar recordings are generated from these unipolar measurements: M1-M2 and M3-M4. The position of the four electrodes of the catheter at recording time is given by the CARTO system and its spatial coordinates were used in order to generate our computational measurements, as shown in Figure 5.5. Simulations took an average of 6 hours to run on a computer including Intel Core i5 CPU and a NVIDIA GT218M(NVS 3100M). Details for each of the five patients are shown in Table 5.2.

Signal Analysis

In addition to a qualitative assessment of signal morphology, a quantitative evaluation of the simulated signals was performed. The following characteristics were extracted from both simulated EGMs and the clinical CARTO recordings:

-Signal Range: Obtained as difference between the maximum and minimum amplitude signal values.

-Number of Inflection Points: Number of points where the signal changes concavity.

-Signal Energy: Calculated using Teager's energy operator [Nguyen 2010].

-Dominant Frequency: Obtained by a fast Fourier transform.

-Mean Slope: Computed as the mean absolute value of dV /dt.

-Fractionation Index: Number of deflections with an amplitude >0.2mV from the signal baseline.

-Minimum Signal Value.

-Maximum Signal Value.

The distributions of the values of each of these characteristics in the LAVA and non-LAVA group were assessed using the non-parametric Kolmogorov-Smirnov (KS) test. A value of p < 0.05 was considered as statistically significant.

Results and Discussion

Catheter measurements were simulated using the image-driven personalized heart model at sinus rhythm. The resulting EGMs were qualitatively compared to their clinical counterparts and features between non-LAVA and LAVA groups were assessed.

CARTO Signal Characterization

Results of the KS test on the signal characteristics of the clinical CARTO recordings are shown in Table 5.4 at the end of this chapter. With the exception of the dominant frequency in Patients 1 and 4, a statistically significant difference between the signal characteristics of non-LAVA and LAVA distributions was found. As expected, LAVA signals present a higher fractionation index and a greater mean slope than their non-LAVA counterpart. The number of inflection points and signal peaks were also higher. A similar difference is expected to be found in the simulated signals. For the non-LAVA simulated signals, we observe a normal ventricular depolarization in the form of a steep negative slope, without any fragmentation inside the EGM signal. The unipolar EGM does not show fragmentation at all. For the LAVA signals, fractionation is seen as a small component with max negative slope. This fractionation could be seen anywhere during the QRS complex of the bipolar EGM. Interestingly, we could also easily detect this fractionation on the unipolar EGM signals. The difference with clinical LAVA is meanly the amount of fractionation and the fact that we find larger LAVA in the real world than in the simulated EGMs.

Simulated EGM Characterization

Results of the KS test for the simulated signals are included in Table 5.3 and a graphic representation is shown in 5.11. The distributions of mean slope, number of inflection points and maximum signal value differed significantly among simulated non-LAVA and LAVA signal populations, as can be seen by the low p value obtained after the KS test. Figures 5.12 and 5.13 show the boxplots of the distributions for two of the features which had the best separation in simulated signals.

Conclusion

We presented the use of a personalized DE-MRI image-based model for the simulation of intracardiac electrograms with LAVA-like characteristics. The phenomenological Mitchell-Schaeffer model was used with literature parameters for the healthy myocardium and appropriate modifications were done in the scar and grey zone to account for the physiological characteristics of these tissues. A dipole approach to intracardiac EGM simulation was also implemented.

The aim of the current work was to explore the feasibility of using in silico models constructed with patient-specific non-invasive data to simulate distinguishable LAVA and non-LAVA intracardiac EGMs. Their resemblance to clinical CARTO signals was qualitatively assessed. Furthermore, feature extraction was performed on both the simulated and clinical signals. Characteristics such as the mean slope and number of inflection points presented significantly different distributions

Introduction

Radio-frequency ablation (RFA) target identification for the treatment of lifethreatening cardiac arrhythmias remains an open challenge. The complete elimination through ablation of local abnormal ventricular activities (LAVA) has been associated with a superior clinical outcome [START_REF] Jaïs | [END_REF]]. The generation of this abnormal electrical activities have been related to regions with the presence of surviving myocardial fibers within scar tissue. Scar topology and tissue heterogeneity can be assessed using non-invasive imaging measures such as delayed-enhanced magnetic resonance imaging (DE-MRI). In previous works [Cabrera-Lozoya 2014], we demonstrated the feasibility in the identification of RFA targets using solely DE-MR image-based features in a machine learning framework while integrating the influence of the inherent error sources of the complex multi-modal data in the training set. Furthermore, we presented in Chapter 5 that the simulation from DE-MR images of distinguishable healthy and LAVA-like intracardiac electrograms.

The aim of this work is to develop an enhanced learning scheme using real and synthetic data to guide RFA therapy by improving the identification of ablation targets, which are defined by the presence of LAVA. First, the performance of the classification is assessed when using solely real or synthetic data. Then, three different approaches merging both data sources are evaluated. The highest performance was obtained when using a fused-learning algorithm. The results derived from this study open up possibilities for non-invasive cardiac arrhythmia ablation planning.

Clinical Data

Five patients referred for cardiac ablation for post-infarction in ventricular tachycardia were included in this study. The patients underwent cardiac MRI prior to high-density EP contact mapping of the endocardium.

Electrophysiological Data. The CARTO mapping system (Biosense Webster) enables the 3D localization of the catheter tip and provides the distribution of EP signals on cardiac surfaces. Contact mapping was achieved in sinus rhythm on the endocardium (trans-septal approach) with a multi-spline catheter (PentaRay, Biosense Webster). Signals were categorised as normal or LAVA by an experienced electrophysiologist.

Imaging Data. The scar tissue was imaged on a 1.5 Tesla clinical device (Avanto, Siemens Medical Systems) 15 minutes after the injection of a gadolinium contrast agent. A whole heart image was acquired using an inversion-recovery prepared, ECG-gated, respiratory-navigated, 3D gradient-echo pulse sequence with fat-saturation (1.25×1.25×2.5mm 3 ). The myocardium was manually segmented on reformatted images of isotropic voxel size (0.625mm 3 ). Abnormal myocardium (dense scar and grey zone areas) was segmented using adaptive thresholding of the histogram, with a cut-off at 35% of maximal signal intensity. Segmentations were reviewed by an experienced radiologist, with the option of manual correction.

Methods

Learning from Image Features

Image based features included intensity and texture features. Intensity features included minimal, maximal, mean and standard deviation intensity values of voxels contained within the catheter's sensing range. Myocardium thickness, scar transmurality and standard deviation over the average intensity in the region were also considered. Haralick features [Tesar 2007] of the myocardial tissue neighboring the site of the electro-anatomical point (EAP) projection onto the endocardium were used to describe the tissue texture within the region of interest. Concatenation of the intensity and texture features yielded a final image-based feature vector of 475 dimensions which was used for classification.

We used a random forest [Criminisi 2011a] classification framework with a nested cross-validation scheme [Ruschhaupt 2004] and optimized the classifier for precision performance using the scikit-learn library [Pedregosa 2011]. Furthermore, we assigned a confidence weighting value to the samples during the training phase based on their temporal displacement during EGM recording. This way, image features from EAPs which were less affected by movement were considered more reliable for the forest construction. The methods have been more thoroughly described in Chapter 4.

Learning from Simulated Intracardiac EGM

Imaging data was used to construct a personalized bi-ventricular model, including distributions of scarred and grey zone regions. The GPU implementation of the Mitchell-Schaeffer model in the SOFA public framework was used along with tissuespecific characteristics to simulate cardiac electrophysiology.

A dipole approach was used for unipolar electrogram simulation at locations where clinical recordings were performed, these spatial coordinates were obtained from the CARTO system. Bipolar electrograms were later computed as the difference between the unipolar measurements of two consecutively placed electrodes.

Feature extraction was performed on the distal-most bipolar simulated EGMs. The feature vector was composed of eight elements: signal range, number of inflection points, signal energy, dominant frequency, mean slope, fractionation index and minimum and maximum signal values. Details concerning the simulation of intracardiac electrograms and EGM feature extraction were discussed in Chapter 5.

The same random forests classification framework with a nested cross-validation scheme as the one described for image-based learning was used. Similarly, the classifier was optimized for precision performance. Nonetheless, no sample weighting Chapter 6. RFA Target Prediction: Combining Imaging Data and Biophysical Modeling due to temporal displacement was included as the precise locations of the catheter electrodes are known during simulation.

Fused-Learning Schemes

Three fusion schemes were proposed to combine the MR image-based and simulated EGM-based learning frameworks:

-Maximum Confidence Fusion. Classification results were obtained independently using image-based or simulated EGM-based features and the confidence in their predictions were analyzed. The final label assigned to the EAP was that of the learning scheme with the highest prediction confidence.

-Image-based, Simulation-Enhanced Fusion. This learning fusion scheme gives greater importance to the classifications performed through image-based learning with temporal displacement weighting. Classification results obtained through simulated EGM-based learning were used solely in the instances where the image prediction had a confidence level below 0.75. In this case, the simulated EGM-based prediction was analyzed. If it was in disagreement with that of the image-based classification and had a higher confidence label, the prediction was modified to that given by the simulated EGM-based classifier.

Else, the image-base classification was left intact.

-Feature Set Fusion. A new feature set was created using both the imagebased and the simulated EGM-based features that have been described in the previous section. This combined feature set was fed to the random forest classification framework using temporal displacement weighting.

Evaluation Metrics

The results of the classification algorithms were projected onto an endocardial surface mesh to create confidence maps for potential RF ablation targets. The study in [Ilg 2010] reports an average endocardial area of 3.5cm 2 for RFA lesions. Therefore, the endocardial surface was divided into regions of area 1.7cm 2 , two times smaller than the average RFA lesions. Figure 6.1 shows a sample endocardial surface mesh partition.

EAP were projected to their closest endocardial surface region along with their prediction and confidence results. Also, a region was considered LAVA if at least one of the EAP projected onto it was labeled as such. The justification behind this is given by the physical constraints of RFA: the ablation of a LAVA site will affect neighboring tissue within the extension of the RFA lesion size, even if they are considered non-LAVA inducing.

After projection of the classification results onto the regions in the surface mesh and its comparison with the ground truth, a confusion matrix was generated and the following statistics were computed: -Precision, or positive predictive value, is intuitively the ability of the classifier not to label as positive a sample that is negative and is computed as P P V = T P/(T P + F P )

where T P, T N, F P, F N correspond to the true positive, true negative, false positive and false negative values of the confusion matrix, respectively.

Results and Discussion

This section presents the results from the previously described learning schemes and ablation target maps from Patient 1 are used to illustrate the algorithm's performance. Further statistics on the classification results are shown in Table 6.1 and Figure 6.7.

Results from Image-based Learning

Classification maps for Patient 1 with their respective confidence values for the LAVA region predictions obtained through MR image-based features are shown in Figure 6.2 for illustrative purposes. As can be seen from Table 6.1, the algorithm has an overall accuracy of 92.9% across the five patients, with a LAVA-specificity of 99.4% and a LAVA-precision of 96.4%. This means that when a region is classified as LAVA, in general, the algorithm is quite confident. Nonetheless, there is a severe compromise on the LAVA-sensitivity, which has an overall score of 33.2%. Classification maps for Patient 1 with their respective confidence values for the LAVA region predictions obtained through the classification of simulated intracardiac electrograms are shown in Figure 6.3. Although the results for this prediction scheme are slightly lower in terms of LAVA-specificity (97.8%) and LAVA-precision (89.4%) than those obtained through MR image-based features, they remain high. The overall accuracy of the algorithm across the five patients is of 94.4%, but more importantly, the LAVA-sensitivity presents a considerable increase, with a value of 73.2%.

Results from Simulated EGM-based Learning

Results from Maximum Confidence Fusion

Results from Image-based, Simulation-Enhanced Fusion

Results for the confidence maps of Patient 1 are presented in Figure 6.5. Although the accuracy achieved is the same as with the maximum confidence fusion scheme and there are slightly higher scores for LAVA-specificity (99.8%) and LAVA-precision (98.8%), there is a drop in LAVA-sensitivity (57.4%).

Results from Feature Set Fusion

The graphical results for this fusion scheme on Patient 1 are shown in Figure 6.6. This method obtained both the top accuracy (97.2%) and LAVA-sensitivity (82.4%) scores while yielding high values for LAVA-specificity (99.2%) and LAVA-precision (95.0%) across the five patients. Therefore, it can be said that this scheme 6.8 shows the mean slope value for these signals. They lie between the 25 th and 75 th percentiles values expected for a LAVA signal, dissipating the suspicions of these signals having non-typical or outlier LAVA signatures.

Classification Results Summary

False Positives Evaluation

In Patient 1, from Figure 6.6 it can be seen that one endocardial region was initially misclassified by the classification algorithm as a false positive and presented a low confidence level (0.65).By analyzing the characteristics of the electroanatomical Figure 6.8: Verifying ground truth point, it was found that the considered sample was one where the catheter presented a high degree of temporal displacement. Its normalized weight was of 0.54, where the weights range from 0.5 for the EAP with the highest degree of movement and 1 for those with most stable positions. So far, sample weighting has been incorporated into our learning scheme in the training phase, but not during the testing phase. This exemplifies the need to consider sample weighting during testing phases of our classification algorithm.

Conclusion

We presented the use of a machine learning framework using real and synthetic data towards the prediction of RFA target identification, defined by the presence of LAVA. The real data was obtained from delayed-enhanced MR imaging while the synthetic was produced using a personalized image-based model for the simulation of intracardiac electrograms. First, the performance of image-based versus simulationbased learning algorithms was compared. Next, several fused-learning frameworks were assessed. The highest performance scores were obtained when both feature sets were merged into a single dataset, the feature set fusion scheme, yielding a mean 97.2% accuracy and 82.4% LAVA-sensitivity scores across five patients. Error regions were characterized by containing EAP with high temporal displacement values. The results derived from this study open up possibilities for non-invasive cardiac arrhythmia ablation planning.

Chapter 7 In this thesis, we presented structurally and functionally personalized models to study the spatial heterogeneity of electrophysiological tissue properties and their role in arrhythmogenesis. Then, we demonstrated that it was possible to use delayedenhanced MR image features to accurately predict potential cardiac ablation targets, based on the presence of LAVA signals. Additionally, we presented a framework for the simulation of intracardiac electrograms using a personalized image-based model. This led to the generation of EGMs with distinguishable LAVA and non-LAVA class characteristics. In the previous chapter, we developed a combined approach to enhance the prediction of RFA targets, using both personalized simulated EGMs and real MRI patient data. In this chapter, we summarize the contributions of each chapter presented in this thesis and propose future directions of research.

Conclusion and Perspectives

Contributions

VT Inducibility Prediction: A combined modeling and clinical approach

The work in Chapter 3 integrated personalized structural and functional data for the study of spatial tissue heterogeneity of electrophysiological (EP) properties using the Mitchell-Schaeffer simplified biophysical model. We showed that patients with highly heterogeneous EP properties were at higher risk of presenting VT. More precisely, cardiac regions with high heterogeneity of restitution properties, such as apparent conductivity and APD-RC slope, were found at the locations of VT exit points. Additionally, we validated the importance of the use of personalized parameters for an accurate in silico VT stimulation procedure by using literature-based parameter values, which either did not induce VT or induced a different macroscopic VT circuit morphology from that observed in clinics. This work was submitted as [Cabrera-Lozoya 2015c].

Image-based Prediction of Cardiac Ablation Targets

The work in Chapter 4 presented a machine learning framework, using intensity and texture-based local imaging features in the vicinity of myocardial scar and grey zones towards the prediction of RFA target localisation. Additionally, we detailed the uncertainty in the data and explored its impact on the classification results.

We proved the feasibility of classification using solely image-based features on complex multi-modal data using three patients. Furthermore, we demonstrated that results were improved when using uncertainty evaluation to weight the training data. Finally, visual interpretation of our results in the form of confidence maps in an anatomical (endocardial or epicardial) surface was presented which has potential use in a clinical environment. This work was published in [Cabrera-Lozoya 2014] and in [Cabrera-Lozoya 2015d].

Image-based Simulation of LAVA Intracardiac Electrograms

In Chapter 5, the use of a personalized DE-MR image-based model for the simulation of intracardiac electrograms with LAVA-like characteristics was presented. The phenomenological Mitchell-Schaeffer model was used with literature parameters for the healthy myocardium and appropriate modifications were done in the scar and grey zone to account for the physiological characteristics of these tissues. A dipole approach to intracardiac electrograms (EGM) simulation was also implemented.

We succeeded in using in silico models constructed with patient-specific noninvasive data to simulate distinguishable LAVA and non-LAVA intracardiac EGMs. Signal features such as mean slope and number of inflection points presented significant differences between the LAVA and non-LAVA distributions for both the clinical and simulated signals. This work will be submitted as [Cabrera-Lozoya 2015a].

Additionally, a Python-based graphical user interface (GUI) was developed in order to obtain a better understanding of the relationshop between the LAVA electrograms and the involved tissue. With this GUI, it is possible to simultaneously visualize the position of the recording electrodes, the image patch associated to this location and to compare the corresponding clinical and simulated electrograms. A screenshot with superposed labels of the functionality of the GUI windows is shown in Figure 7.1.

RFA Target Prediction: Combining Imaging Data and Biophysical Modeling

In Chapter 6, the MR image-based machine learning framework described in [Cabrera-Lozoya 2014] was combined with the personalized intracardiac EGM simulation scheme in [Cabrera-Lozoya 2015a] to enhance the prediction of RFA target identification. We described how the inclusion of simulation data improved the performance of cardiac ablation target identification obtained from imagebased learning. The results derived from this study open up possibilities for non-invasive cardiac arrhythmia ablation planning. This work will be submitted as [Cabrera-Lozoya 2015b].

Methodological and Clinical Perspectives

Methodological Perspectives

Throughout this thesis, cardiac electrophysiology simulation was performed using the mono-domain Mitchell-Schaeffer [Mitchell 2003] model. Nevertheless, the bidomain model is known to be more accurate in mathematically describing cardiac electrical activity, with a major drawback of being extremely computationally expensive. The work in [START_REF] Keller | [END_REF]] studied the differences in intracardiac signals when using mono and bidomain models and found that bidomain simulations better represent in-vivo data when investigating small-scale changes in intracardiac electrograms. [Coudiere 2014] has developed a computationally efficient and accurate approximation of the bi-domain model that allows for extracellular stimulation and accounts for unequal anisotropy ratios between intra-and extra-cellular media. This opens the possibility of exploring the use of a more refined electrophysiology model in our work while keeping computation times and loads to an acceptable level. Furthermore, in the current state, the EGM simulation implementation could highly benefit from a parallelization setting in a GPU framework, like the one developed in [Talbot 2014] for EP simulations. This would achieve shorter computation times than the ones reported in Chapter 5 and which could be acceptable to use in a clinical environment.

One of the main contributions of this work was the consideration of sample uncertainty during the training phase when electro-anatomical points were weighted according to their degree of temporal movement. Nonetheless, an aspect that was not explored during this work was the role of uncertainty in the test samples. An extension of this work can include the knowledge of uncertainty into the testing phase.

Classification labels of EAP with high temporal displacement can be complemented by additional EGM simulations and image feature computation in locations of the displacement trajectory of the EAP. This would yield a higher number of descriptors that would better characterize the underlying tissue region swept by the recording catheter and could help increase the classification performance by considering the density of classified labels of these supplementary points.

In this work, per-patient analysis were performed for the prediction of cardiac ablation targets. In order to evolve from per-patient analysis into an inter-patient learning framework, a large database should be constructed in order to allow for a varied and diverse training pool. Nonetheless, because clinical data is sometimes hard to obtain, the database size can be increased by the incorporation of synthetic but realistic cases as was done in the work that will appear in [START_REF] Duchateau | Predictions of infarct localization from myocardial deformation[END_REF]]. These samples could include realistic myocardial infarctions with different grey zone to scar core ratios, sizes, shapes and locations and intracardiac electrograms could be computed from them. It is only after considerably increasing the size of our database (with both real and synthetic samples) that we might be able to consider more sophisticated learning algorithms, such as deep learning, which has recently been used in other computer vision applications [START_REF] Taigman | [END_REF]].

The use of new image and signal-based features can be explored. Examples of families of image features that have not been investigated in this work are wavelets or steerable filters. Also, the use of image features with rotation invariant properties would be of great use in an inter-patient study. For the EGM-based counterpart, time-frequency domain features can be incorporated. Other methods have been proposed for the analysis of EGM signals, including those with wavelet decomposition approaches. The study in [Alcaine 2014] proposes the wavelet decomposition on the signal envelope of EGMs for the automatic generation of activation maps. Perhaps the use of signal features in different wavelet scales can help improve the classification of simulated intracardiac electrograms.

Nonetheless, we think that this work would greatly benefit from the use of principal component analysis (PCA) or linear discriminant analysis (LDA) for dimensionality reduction. Retrieving a reduced subset of features without compromising classification performance would render the classifier more robust against overfitting and would also help to speed up the feature extraction phase, particularly the imagebased. A further reduction in computation time can be obtained by computing only the required image features on the fly at each stage of the testing phase.

Clinical Perspectives

A considerable amount of clinical data during this thesis was obtained through the use of a PentaRay catheter, which is not the standard used in clinics. This multispline recording catheter allows for a higher-density mapping in a shorter time lapse. While our framework would not change if a multi-spline catheter is not used, this would considerably increase the clinical procedure time, as the catheter would need to be swept across a larger number of positions in order to obtain a high-density electrical map. In this scenario, previously simulated intracardiac electrograms from image-based patient-specific models could be of great help, by generating maps containing the regions which are more likely to present LAVA electrograms.

In the last years, electrocardiographic imaging systems have been presented as a method for mapping the electrical activity of the heart [START_REF] Rudy | Noninvasive electrocardiographic imaging of arrhythmogenic substrates in humans[END_REF]]. While this technique might be useful when studying large scale activation patterns and it could possibly find applications in the study of ventricular tachycardia, it still lacks the properties to study local phenomena, like LAVA.

Image normalization and image quality assessment are issues that should be managed. Indeed, image quality can play a very important role in the image-feature extraction phase. For example, using training image features from high quality MR images will not perform as expected on a test image with low resolution, and viceversa. Therefore, a means of evaluating image quality should be derived and appropriate classifiers should be developed for each scenario. Furthermore, the usefulness of medical imaging alternatives, such as computed tomography (CT), should also be assessed of our particular clinical problem. In the pre-clinical study described in [Schuleri 2009], the authors compared the characterization of the border zone when using delayed-enhanced multi-detector CT to the more widely used DE-MR imaging techniques. They claimed that DE-MDCT provided a more detailed assessment of the border zone and was less susceptible to partial volume effects than MRI. Therefore, it would be interesting to evaluate the performance of the algorithms proposed in this thesis when using CT imaging as an input.

As a long-term perspective and after further validation of the proposed techniques on a larger patient cohort, we can think of a fully coupled learning and clinical ablation target identification approach as shown in Figure 7.2. The patient will be firstly referred to an MR imaging study from which a patient-specific image-based model could be derived. A preliminary ablation target map would be generated using the developed techniques with an an offline classification stage. Then, during the clinical study, real catheter recordings and the clinician's expertise can be incorporated to complement and update the target map in real time. Because random forests are efficient algorithms that can easily benefit from parallelism, they become ideal candidates for an online learning framework. The simplest approach requires to push newly obtained training data through all the trees to update the leaf distributions, keeping the learned parameters and the forest structure fixed, as has been described in [Criminisi 2011a].

Furthermore, ablated regions can be incorporated in the image-based model and changes in cardiac electrophysiological behavior can be assessed. This would lead to a fully integrated RFA planning scheme for therapy guidance.

In this PhD, we explored how the close integration of mathematics, computer science and medicine can leverage on important progress made in computer simulation and machine learning to improve cardiac intervention planning.

A.1 Introduction

The inherited long-QT syndrome (LQTS) is an arrhythmogenic disease characterized by impaired cardiac repolarization that clinically manifests with QT prolongation, polymorphic ventricular tachycardia (pVT), syncopes, and sudden cardiac death (SCD) [Roden 2008]. A spatially heterogeneous prolongation of repolarization leading to an enhanced dispersion of action potential duration (APD) is considered a major contributor to LQT-related arrhythmias [Antzelevitch 2007]. As a result of their prolonged repolarization, patients with LQTS have globally prolonged cardiac contraction duration with transmural differences [Haugaa 2009], [Haugaa 2010]. However, the spatial relationship between electrical and mechanical (dys)function remains to be elucidated.

Phase contrast magnetic resonance imaging (MRI) has been successfully used to determine regional differences in myocardial contraction and relaxation velocities in human patients [START_REF][END_REF]]. Healthy women, who physiologically have longer repolarization than men, have altered peak velocities and different regional patterns of relaxation than men [Föll 2010], suggesting that it may be possible to determine changes and spatial differences in mechanical function in LQTS using this technique.

Rabbits exhibit pronounced similarities in the expression and function of voltagegated K + channels [Nerbonne 2000] and Ca 2+ cycling proteins [Bers 2002] as humans, suggesting that rabbits can fairly appropriately mimic the human phenotype of cardiac diseases with electrical and/or mechanical impairment. Indeed, several transgenic rabbit models of human diseases such as LQTS and hypertrophic cardiomyopathy have demonstrated pronounced similarities with the human phenotype [Brunner 2008], [Ripplinger 2007]. Transgenic LQT2 rabbits mimic human LQTS with QT prolongation, pVT, SCD, and an increased APD dispersion as major arrhythmogenic mechanism [Brunner 2008]. Moreover, despite some species differences in myocardial motion between humans and small animals, we could recently show similarities between humans and rabbits also in regional contractile behaviour with a similar rotational motion [Jung 2012], indicating that LQT2 rabbits may be a useful tool to assess mechanical function in LQTS.

Using in vivo phase contrast MRI, ex vivo monophasic action potential (MAP) measurements and machine learning approaches, we aimed at elucidating the spatial relationship between regional electrical and mechanical cardiac function and their spatial dispersion in transgenic and drug-induced LQT2 rabbits.

A.2 Methods

A more detailed method description can be found in the online supplement in http: //dx.doi.org/10.1016/j.hrthm.2013.07.038..

A.2.1 Rabbits

Adult female transgenic LQT2 (HERG-G628S, 4.1±0.6 months, 2.4±0.3 kg) and wildtype littermate control (LMC) rabbits (without drugs, 6.3±1.0 months, 3.8±0.6 kg; with E4031, 4.0±0.2 months, 2.9±0.3 kg) were subjected to MRI followed by MAP measurements.

A.2.2 Phase contrast MRI

To assess myocardial velocities, transgenic LQT2 (n=11) and LMC rabbits (n=9 without drugs, n=10 with E4031, bolus 10µg/kg and infusion 1µg/kg/min [Michael 2007]) were subjected to phase contrast MRI at 1.5 T (Avanto, Siemens, Germany) [Jung 2012].

Animals were anesthetized with s-ketamine/xylazine (12.5/3.75 mg/kg IM, followed by 1-2.5 ml/kg/h IV), which does not affect cardiac repolarization [Odening 2008], and positioned in a 12-channel head coil.

Experiments were performed with a black blood prepared gradient echo sequence with prospective ECG-gating and high temporal (7.6 ms) and spatial resolution (1.0×1.2×4 mm) [Jung 2012]. For data post-processing, customized Matlab software was used. For regional analysis, LV was partitioned into 16 segments in base, mid, and apex (AHA model [Cerqueira 2002]). Mean radial (Vr) and long-axis (Vz) systolic and diastolic peak and time to peak (TTP) velocities were determined. Standard deviation (SD) of TTP velocities from 16 segments was calculated as measure for mechanical dispersion. Heart rate corrected diastolic TTP were computed as percentage of RR.

A.2.3 MAP measurements in Langendorff-perfused rabbit hearts

To correlate electrical and mechanical function, MAPs were acquired in the same LQT2 (n=10) and LMC rabbits (n=9 without drugs, n=10 with E4031, 0.1µM [D'Alonzo 1999]). Rabbits were anesthetized with s-ketamine/xylazine (12.5/3.75 mg/kg) IM, and received 1000 IU heparin IV. Following euthanasia with thiopentalsodium (40 mg/kg) IV, hearts were rapidly excised and mounted on Langendorffperfusion set-up (IH5, Hugo Sachs Electronic-Harvard Apparatus, Hugstetten, Germany). Hearts were retrogradely perfused with modified Krebs-Henseleit solution [Bentzen 2011]. ECG, MAP, coronary flow and pressures were continuously recorded with Isoheart software (Hugo Sachs Electronic, Version 1.1.1.218(32)).

After mechanical atrioventricular node ablation, hearts were stimulated with 2, 3, and 4 Hzto obtain APD at heart rates comparable to MRI (130-189 beats/min). Four MAP electrodes were repetitively positioned on all different LV segmentsexcept for the septal segments (Figure A.2). APDs at 75% of repolarization (AP D 75 ) and standard deviations (SD) of AP D 75 within all LV segments were calculated as measure for electrical dispersion. Transgenic and Drug-induced LQT2 Rabbits

A.2.4 Machine learning for image-based classification

Multivariate data analyses were performed using random forest machine learning algorithms to explore the possibility to differentiate between LQT2, E4031, and LMC based on MRI. Image-based features included peak velocities, time-to-peak velocities, regional velocity statistics and correlations to normal velocity curves (over 700 features/rabbit). Feature space dimension reduction was performed through chisquared test. Only top-relevant features were considered. Cross-validation was used to obtain an estimate of classifier accuracy [START_REF] Menze | [END_REF]].

Since heart-rate corrected QT (QTc) duration [Sauer 2007] and QT dispersion [Priori 1994] are known arrhythmogenic risk factors in LQTS, and since we observed correlations between APD and diastolic dysfunction in LQT2 rabbits, we secondly tested whether MRI can provide sufficient information to discern differences in the extent of APD prolongation using random forest algorithms. LQT2 rabbits were divided in two groups: 1) rabbits with "very long" APD and "pronounced" dispersion (longest APD > 155 ms, SD_APD > 12 ms) and 2) rabbits with "long" APD and "moderate" dispersion (longest APD > 140 ms and < 155 ms, SD_APD < 12 ms) (Table 1 suppl.).

A.2.5 Statistical Analysis

For normally distributed values, we used unpaired Student's t-test for comparisons and Pearson's test for correlations. For values not normally distributed, we used Mann-Whitney test for comparisons and Spearman's test for correlations. Correlations were visualized using linear least square fit calculated based on all LQT2 rabbits' data. Analyses were performed with Prism 4.03 (Graphpad, San Diego, USA).

A.3 Results

A.3.1 Dispersion of Cardiac Repolarization in Transgenic and

Drug-induced LQT2 Rabbits During faster stimulation, APD shortened in all groups (LQT2, p<0.0001, E4031, p=0.002, LMC, p=0.002). Despite this rate-dependent APD shortening, however, APD remained significantly longer in LQT2 and E4031 than in LMC at 3 Hz (3Hz, base, ms, LQT2, 132±8.5, p=0.02, E4031, 121±3.6, p=0.04, LMC, 113±1; 4 Hz, LQT2, 105±9.7, E4031, 97±7.6, LMC, 89±0.5). Moreover, APD dispersion remained more pronounced in LQT2 and E4031 than in LMC at 3 and 4 Hz (Figure A.1D).

A.3.2 Mechanical Function in Transgenic and Drug-induced LQT2 Rabbits

Global peak systolic radial velocities (Vr) were significantly reduced in LQT2 and E4031 compared to LMC in LV base (LQT2, p=0.01; E4031, p=0.0004) and mid (LQT2, p=0.002; E4031, p=0.0009), and tended to be reduced in apex (LQT2, p=0.07) -due to regionally reduced velocities in anterior and anteroseptal segments ( Global and regional peak diastolic Vr were significantly reduced in LQT2 and E4031 in LV base (LQT2, p=0.002; E4031, p<0.0001), mid (LQT2, p=0.002; E4031, p<0.0001), and apex (LQT2, p=0.01; E4031, p=0.006; Figure A.3A, B) -indicating an impaired diastolic function in LQT2. Global peak diastolic Vz did not differ between LQT2 and LMC, but were reduced in E4031 (base, p<0.0001; mid, p=0.0002; apex, p=0.02). Regional diastolic Vz were reduced in LQT2 in anterior segments, while in E4031, diastolic Vz were reduced in inferior and lateral regions (Figure A.3D).

Global and segmental time-to-diastolic peak Vr (TTP_dia_Vr) were significantly longer in LQT2 and E4031 than in LMC (LQT2, base, p=0.0002; mid, p=0.04; E4031, base, p<0.0001; mid, p<0.0001, Figure A.4A) -indicating a delayed relaxation due to prolonged contraction duration. After heart-rate correction, differences between LQT2 and LMC were even more pronounced (TTP_dia_Vr_base%,LQT2,65.2%±0.8,LMC,58.8%±0.6,p<0.0001;Figure A.4B). In contrast, global and segmental heart-rate corrected TTP_dia_Vz were similar in all groups.

A.3.3 Correlation of electrical and mechanical dysfunction

In LQT2 rabbits, APD correlated moderately with regional reduction of diastolic Vr and Vz in LV base (Vr_dia,CC 0.38,p=0.02;Vz_dia,CC 0.47,p=0.002;Figure A.3C). Moreover, SD of APD and peak diastolic Vras markers of electrical/mechanical dispersioncorrelated strongly in LV mid (CC 0.75, p=0.03). In addition, TTP_dia_Vr -as marker of contraction durationcorrelated moderately with APD in LV base (TTP_dia_Vr,CC 0.47,p=0.001;TTP_dia_Vr%,CC 0.43,p=0.001,Figure A.4D). As for systolic function, APD correlated moderately with systolic Vz in LV mid (CC 0.37, p=0.03). In E4031-treated rabbits, APD correlated with diastolic Vr (CC 0.42, p=0.04), while in LMC, no correlations were observed.

A.3.4 Machine learning approaches for MRI-based classification

Based on the observed differences in mechanical function, we first explored the possibility to differentiate between groups solely based on "mechanical" phase contrast MRI. Using random forest algorithms with only most relevant features (Table 2 suppl.), we could differentiate between LQT2 and LMC with an accuracy of 95 %. All LQT2 rabbits were correctly classified, while one LMC with relatively low Vr_dia was misclassified (Table 3 suppl.). Most relevant features were regional diastolic Vr, TTP_dia_Vr, and their SD (Table 2 suppl.), while systolic features were only of minor importance. As indicated in the cumulative region relevance plot, important features were mainly located in base (Figure A.5A). Moreover, we could differentiate between E4031 and LMC with an accuracy of 100% (Table 3 suppl.). Most of the relevant features for classification were regional diastolic Vr and Vz, TTP_dia_Vr, and SD of Vr_dia, TTP_dia_Vr, and TTP_dia_Vz (Table 2 suppl.), mainly located in base inferior and mid regions (Figure A.5B).

Based on the correlations between electrical and mechanical features in LQT2, we secondly tested whether MRI could provide sufficient information to discern even slight differences in the extent of APD prolongation (LQT2 "very long" APD, 164±7.6 ms; LQT2 "long" APD, 144±2.8 ms, p<0.001, Table 1 suppl.). Using only the 17 most relevant features (Table 2 suppl.), it was possible to differentiate between"very long" vs. a "long" APD groups with a classification accuracy of 100 % (Table 3 suppl.). Most relevant features were diastolic Vr and Vz and TTP_dia-TTP_sys (Table 2 suppl

A.4 Discussion

Using phase contrast MRI, MAP, and multivariate random forests machine learning algorithms, this study provides a comprehensive assessment of regional dispersion of repolarization and mechanical function in LQT2 rabbits. We revealed that a genetically-induced or I Kr -blocker-induced prolongation of APD led to a globally and regionally impaired diastolic function and, importantly, that regional differences in APD correlated with regional differences in diastolic dysfunction. These observations indicate that inherited and drug-induced LQTS are not purely electrical but rather electromechanical disorders. Moreover, we demonstrated the possibility to a) differentiate between LQT2 and normal rabbits, and b) discern differences in the extent of APD prolongation in transgenic LQT2 rabbits solely based on phase contrast MRI, indicating a potential future use of MRI for risk stratification in LQTS.

A.4.1 Dispersion of Repolarization in LQT2 Rabbits

A spatially heterogeneous prolongation of repolarization leading to an enhanced APD dispersion is considered a major contributor to LQT-related arrhythmias [Antzelevitch 2007]. LQTS patients have a more pronounced QT and Tpeakend dispersion indirectly reflecting regional and transmural APD heterogeneities [Priori 1994], [Lubinski 1998]. Spatial and transmural heterogeneities in I Kr , I Ks , and I Ca,L underlie these pronounced APD heterogeneities in LQTS [Cheng 1999], [Pham 2002]. To directly assess dispersion of repolarization in LQTS patients, invasive electrophysiological studies (EPS) are necessary and this assessment is limited to regions easily accessible during EPS [Bhandari 1985]. Taking advantage of transgenic LQT2 rabbits, we have previously demonstrated a pronounced dispersion of repolarization on RV endocardium [Odening 2010] and LV anterior surface [Brunner 2008] using in vivo EPS and ex vivo optical mapping. We here complemented these data by systematically assessing APD in anterior, anterolateral, inferolateral, and inferior regions of LV base, mid, and apex and demonstrated a more pronounced APD dispersion within each level and global LV. Moreover, we demonstrated that this increased APD dispersion was not altered by faster stimulation despite a rate-dependent APD shortening.

A.4.2 Dispersion of Contractile and Diastolic Function in LQT2 Rabbits

For the first time, we systematically assessed global and regional systolic and diastolic function in LQTS and demonstrated that APD prolongation a) was associated with delayed cardiac relaxation and b) had an impact on radial systolic velocities andparticularlyon radial and longitudinal diastolic velocities, indicating an impaired mechanical function in LQT2.

Heart size (LV volumes, ventricular mass, and wall thickness) may impact on systolicand less pronouncedly on diastolicvelocities [Föll 2010]. The observed pronounced differences in systolic and diastolic velocities in LQT2 and E4031 rabbits compared to LMC, however, cannot (solely) be explained by differences in LV volumes since E4031-treated rabbits were bigger than LQT2 but had similarly reduced diastolic velocities. Furthermore, differences in heart size cannot explain regional correlations between APD and diastolic velocities in transgenic and drug-induced LQT2 rabbits as heart size did not vary much within these groups.

Since electrical and mechanical cardiac functions are coupled and influence each other [ter Keurs 2011], LQTS patients have not only prolonged repolarization but also have similarly altered ventricular contraction patterns with slow contraction in the late wall thickening phase [Nador 1991], longer isovolumetric relaxation times [Savoye 2003], and globally prolonged contraction duration with transmural differences [Haugaa 2009], [Haugaa 2010] as demonstrated by M-mode and tissue Doppler echocardiography. Moreover, in a case report of an infant with a particularly long QT, the onset of pVTs was preceded by a globally impaired systolic and diastolic function [Vyas 2008]. Two small studies demonstrated that mechanical abnormalities in LQTS patients were associated with arrhythmogenic events [Nador 1991]: the alteration of myocardial velocities was more sensitive and more specific for the prediction of cardiac events than QTc duration [Haugaa 2009], indicating a possible prognostic value of regional velocity imaging in LQTS patients.

The tissue Doppler imaging technique, however, is limited by low reproducibil-Transgenic and Drug-induced LQT2 Rabbits ity, angle dependency of measurements, and does not permit the evaluation of all myocardial segments and all velocity directions [START_REF] Thibault | [END_REF]]. Phase contrast MRI overcomes these limitations and allows assessing all velocities in all segments [START_REF][END_REF]]. Using this technique, we revealed that not only global peak velocities are reduced and global time-to-diastolic peak velocities are prolonged in LQT2, but rather that reduction of diastolic velocities and prolongation of contraction duration are regionally heterogeneous, indicating a mechanical dispersion. Interestingly, in transgenic LQT2 rabbits, radial velocities were particularly reduced and radial time-to-diastolic peak were particularly prolonged. In line with our findings, in LQTS patients a very slow late contraction was observed in radial direction [Nador 1991]. Moreover, dispersion parameters and contraction duration in radial velocities demonstrated the highest correlations to APD. These observations stress the particular importance of evaluating regional radial motion in LQTS, which cannot be assessed in all LV segments using standard echocardiography.

A.4.3 Correlation of electrical and mechanical dysfunction and MRI-based differentiation of electrical phenotype

Although it has previously been shown that LQTS patients have a pronounced dispersion of both, cardiac repolarization and (transmural) contraction duration [Haugaa 2010], thus far, no study has investigated a potential spatial relationship between electrical and mechanical dispersion in LQTS. Taking advantage of transgenic LQT2 rabbits [Brunner 2008], in which mechanical and electrical function can be assessed systematically in similar segments, we revealed for the first time a spatial correlation between APD and the extent of diastolic dysfunction, contraction duration, and mechanical dispersion, indicating a close coupling of spatially heterogeneous electrical and mechanical dysfunction in LQTS. Due to electro-mechanical coupling in cardiomyocytes, one might expect an even stronger statistical correlation than the one observed in our study. However, regional differences in fiber orientation and heart geometry at different segments may additionally affect amplitude and timing of radial, longitudinal and rotational velocities thus modulating the direct impact of APD on relaxation characteristics. Nevertheless, this observed association between electrical and mechanical changes in LQT2 rabbits was pronounced enough to discern even moderate differences in the extent of APD prolongation solely based on phase contrast MRI. How this altered electro-mechanical function in LQTS may impact on arrhythmogenesise.g. via reverse excitation-contraction coupling [ter Keurs 2011], via activation of stretch-sensible ion channels [Wang 2009b], or via reduced coronary perfusion due to prolonged cardiac contraction and impaired diastolic relaxation remains to be elucidated. APD prolongation results in increased Ca 2+ influx via I Ca,L and increased Ca 2+ -induced Ca 2+ -release from the sarcoplasmic reticulum, thus increasing cytoplasmic Ca 2+ concentration [Ca 2+ ] i [Bassani 2006]. Since physical properties of (resting) sarcomeres are related to [Ca 2+ ] i [Stuyvers 2000], an increased [Ca 2+ ] i may result in an impaired diastolic relaxation. Moreover, the in-creased [Ca 2+ ] i should additionally increase systolic function31. However, our data suggest that systolic function in LQT2 hearts is decreased, indicative of a complex remodeling of Ca 2+ handling systems like an adaptive decrease in I Ca,L density observed in transgenic LQT2 rabbits [START_REF] Liu | [END_REF]]. Further studies are necessary to fully characterize molecular mechanisms responsible for adaptive changes contributing to systolic and diastolic dysfunction in LQT2.

A.4.4 Clinical implications

To date, LQTS is mainly considered an "electrical" disorder, in which QTc duration [Sauer 2007] and QT dispersion [Priori 1994] are known arrhythmogenic risk predictors. The fact that a) regional differences in APD correlate with regional differences in diastolic dysfunction and that b) it is possible to discern differences in the extent of APD prolongation using phase contrast MRI might open new approaches for risk stratification and treatment decisions of LQTS patients. In the future, it might be possible to use phase contrast MRI to assess the arrhythmogenic risk in LQTS patients non-invasively by deducing electrical dispersion from segmental mechanical function and dispersion. Transgenic LQT2 rabbits overexpress the dominant-negative loss-of-function mutation HERG-G628S that leads to a complete loss of I Kr . Many missense mutations in HERG in LQT2 patients, however, lead to a substantial decrease rather than a complete loss of I Kr . The findings of pronounced diastolic dysfunction in transgenic LQT2 rabbits thus likely recapitulate findings in human patients with a (complete) loss of functional I Kr . The fact that I Kr -blocker E4031, which only reduces I Kr , also induces diastolic abnormalities, however, indicates that LQT2 patients with a substantial decrease of I Kr may have similar diastolic dysfunction. To validate the potential use of phase contrast MRI for risk stratification, further studies investigating its potential to predict the development of pVTs and clinical pilot studies investigating the extent of mechanical (dys)function and its spatial dispersion in LQTS patients with different mutations, different QTc, and different arrhythmia incidences are certainly warranted.

A.4.5 Limitations

Regional mechanical and electrical data were acquired in different experimental settings, e.g. in vivo MRI and ex vivo MAP in Langendorff-perfused hearts. Ideally, both should be assessed simultaneously in the same in vivo setting to perfectly correlate regional mechanical and electrical function. However, to date, detailed segmental MAP measurements are not feasible in an in vivo setting and it is still very challenging to obtain anylet alone detailedregional electrical data during MRI. The "diastolic dysfunction", e.g. abnormal cardiac relaxation, observed in LQT2 rabbits, must be differentiated from "diastolic heart failure". Clinically, none of the LQT2 rabbits showed any signs of diastolic heart failuresimilarly as in LQTS patients, in whom no clinical signs of heart failure were reported, even in older patients or in patients with particularly prolonged QT [START_REF] Goldenberg | [END_REF]]. Moreover, so far, the muscle cells. Since [START_REF] Hodgkin | [END_REF]], several approaches have been developed to describe the electrophysiology, from biophysical models that simulate ionic interaction [Noble 1962], to phenomenological models based on reaction-diffusion equations [FitzHugh 1961] [Aliev 1996] and recent adaptations of those [Clayton 2011[START_REF] Pashaei | [END_REF]] [Relan 2011a]. They differ in the scale they use and their complexity in terms of computation times and number of parameters. These models have to be personalized from endocardial mappings and electrocardiograms [?] before serving as input to the mechanical simulation.

Finally, the electromechanical coupling itself has a high impact on the resulting simulation. For the past 20 years, the creation of such models has been the focus of many research teams and various models have been proposed [START_REF] Humphrey | [END_REF], Hunter 1997[START_REF] Nash | Mechanics and material properties of the heart using an anatomically accurate mathematical model[END_REF], Bestel 2001, Sachse 2004]. In this paper, the Bestel-Clément-Sorine cardiac model [Bestel 2001], further improved by [Chapelle 2012b] is used. It shows a good compromise between complexity and accuracy [Marchesseau 2012b], and a good predictive power for CRT under different pacing conditions [START_REF] Sermesant | [END_REF]]. The personalization of such mechanical models has been tackled by different approaches. Variational assimilation methods [Delingette 2012[START_REF] Sainte-Marie | Modeling and estimation of the cardiac electromechanical activity[END_REF] were used to locally estimate the active parameters (contractility, contraction and relaxation rates), Sequential Quadratic Programming [Wang 2009a] to estimate the passive material stiffness, while Reduced Order Unscented Kalman Filtering was aimed at estimating the contractility parameters [Chabiniok 2011, Moireau 2011].

These methods have already led to promising results to estimate some model parameters from heart motion and therefore open the possibility to help cardiologists in planning their therapy. In this paper we chose to tackle this issue in a different way. We propose to estimate the regional contractilities of the left ventricle based on measured Regional Volumes. Using regional volumes has the advantage of filtering registration errors by only considering displacements orthogonal to the surface thus accounting for the aperture problem. Moreover, automatic estimation of regional volumes is already available in commercial software which may help to streamline the personalization process. We apply this strategy on a database consisting of eight healthy subjects and three heart failure patients. To this end, we first calibrate the mechanical model using the Unscented Transform method proposed in [Marchesseau 2012a] from the left volume curve that we extract from the motion registered from cine MRI data. Then we apply the Reduced Order Unscented Kalman Filtering [Moireau 2011] -using the Verdandi library 1 -on the regional volumes measured on the LV endocardium. This strategy was successfully tested synthetically for contractility estimation and scar tissue detection, and on real cases. Therefore, the use of regional volumes to personalize the regional contractilities proved to be efficient and robust and testing the full pipeline on 8 healthy and 3 pathological cases enabled to draw preliminary conclusions on parameter specificity. In this re-using Reduced Order UKF from Regional Volumes spect and to the best of our knowledge, this study exhibits the largest cohort of healthy and pathological cases for personalized cardiac models to date.

B.3 Materials and Methods

B.3.1 Data Acquisition

The proposed study was performed on 3 pathological and 8 healthy cases. Data were acquired either at INSERM (Institut National de la Santé et de la Recherche Médicale, Université de Rennes 1, France), or at the Division of Imaging Sciences & Biomedical Engineering at King's College London, UK, as part of studies that were ethically approved. Written informed consent was obtained from all human subjects and the studies were performed with all appropriate institutional human subjects approvals. All datasets consist of sequences of cine SSFP MRI with a spatial resolution of 1.5x1.5x7 mm 3 for the healthy subjects and 1.25x1.25x10.8 mm 3 for the pathological cases. Moreover, an intra-operative electrophysiological study (EnSite) was performed on all pathological cases.

B.3.1.1 Real cases study

All heart failure patients were selected for Cardiac Resynchronization Therapy based on criteria from the New York Heart Association (NYHA) class III-IV. Patients had a small ejection fraction (EF) < 35% and a prolonged QRS on electrocardiogram > 120ms.

This study also includes extensive multi-modality imaging of volunteers from which 8 healthy cases were used. Volunteers were aged 28 5 years, without clinical history of cardiac diseases. This database was made available to the research community for 1st Cardiac Motion Analysis (cMAC) challenge held at STACOM'2011, see [Tobon-Gomez 2013] for details regarding the data acquisition of this study.

B.3.1.2 Clinical Data Pre-Processing

Three different steps are needed before any mechanical personalization can be performed: extraction of the myocardium geometry, estimation of the patient's cardiac motion and personalization of the electrophysiological propagation.

Geometry Personalization

We first extract a biventricular anatomy from SSFP MRI datasets using the segmentation algorithm implemented in [Tobon-Gomez 2012]. This semi-automatic segmentation technique is based on a deformable model and allows manual deformations if final corrections are needed. We generate a volumetric binary image from this biventricular surface (see Fig. B.1) after interpolation and resampling on the zaxis for a higher resolution mesh. The Marching Cubes algorithm [START_REF] Lorensen | [END_REF]] is then used to extract a polygonal surface. ReMesh [Attene 2006] enables to smooth and ensure a manifold triangular mesh while Tetgen2 and Netgen3 generate the tetrahedral meshes (see Fig. B.1) with approximately 14 000 nodes. The sensitivity of the mesh resolution on the resulting motion was studied and enabled to conclude that this resolution is a good trade-off between computation times and accuracy with respect to the available data.

In addition, fibre directions were created synthetically as proposed by Streeter [Streeter 1979]. It is derived from both histological and DTI studies [Streeter 1979, Hsu 2001] and it has been extensively validated for cardiac simulation [Bayer 2012]. The elevation angle was varied (with respect to the short axis plane) across the myocardium wall from -70°on the epicardium to 0°at mid-wall to +70°on the endocardium. This choice influences the global motion of the myocardium and was therefore set after a trial and error approach. In the myocardium, the angle was interpolated using the distance from the selected node to the closest surface. Alternatively, a way of creating the fibres could be to map a statistical atlas of human cardiac fibres such as [START_REF] Lombaert | Statistical Analysis of the Human Cardiac Fiber Architecture from DT-MRI[END_REF]] onto the patient geometry. This would allow to include real patient data better than with predefined synthetic angles. However, statistical fibre models have not yet been validated for cardiac simulation.

Motion Estimation

Motion estimation was performed from cine MRI data using an automatic registration algorithm: the Time Diffeomorphic Free Form Deformation (TDFFD) algorithm, which has been tested on 3DUS datasets [De Craene 2012a] and 3DTAG MRI datasets [De Craene 2012b]. The TDFFD algorithm optimizes a 4D velocity field parametrized by B-Spline spatiotemporal kernels.

With the transformation obtained from the registration, the static volumetric mesh was deformed and local or regional indices were computed. Due to its temporal consistency, this technique computes smooth volumetric meshes which is highly desirable for our application.

Resulting meshes were compared to the images in and Fig. B.2 that information on the basal area is very incomplete. Since we cannot estimate thoroughly the real movement in this area, we ignored the observations extracted from the registration related to this area when personalizing the model. This results in smaller regional volume zones (see Section B.3.2.3) but this does not alter the convergence of the algorithm. Moreover this prevents unrealistic parameters to be estimated without influencing the motion drastically.

Electrophysiology Personalization A large number of electrophysiological models for cardiac cell action potential have been developed, ranging from models at the cellular to those at the sub-cellular scale. They can be divided into three categories:

• Biophysical models. Are the most complete because they model ionic currents.

Nonetheless, their high complexity, computational cost and lack of observability of their parameters make them less suitable for parameter estimation from clinical data.

• Phenomenological models. Are based on PDEs and present an intermediate level of complexity and computational cost.

• Generic models. Are the least complex set of models and represent simplified action potentials or depolarization times. Eikonal models belong to this category and model the propagation of the time at which a given point is depolarized.

To simulate the electrophysiological pattern of activity, an Eikonal equation was used. The model was then solved for the depolarization time T d at each point of The local electrical conduction velocity is v and D is the anisotropic conductivity tensor defined as

D = (1 -r)f ⊗ f + rI (B.2)
where f is the fibre orientation, r the anisotropy ratio and I the identity matrix.

The solution of this electrophysiological model was found using a Multi-Front Fast Marching Method [START_REF] Sermesant | [END_REF]]. For healthy cases no subject specific electrophysiological data were acquired and therefore standard values were assumed (conduction velocity of 90cm/s and anisotropic ratio of 0.1). The personalization of this model requires to specify the onset of the electrical propagation (manually defined as areas corresponding to the extremities of the Purkinje network, on both endocardiums close to the apex).

For pathological cases, intracardiac non contact mappings allowed to estimate the time and localization of the onset of depolarization as well as electrical conductivities (see [Relan 2011a]). For the personalization, the Mitchell-Schaeffer [Mitchell 2003] biophysical model is used. It simulates the evolution of the transmembrane potential as a function of simplified inward and outward currents. An example of personalized electrophysiology is given in

B.3.2 Mechanical Personalization Strategy

In this paper, a SOFA4 implementation of the Bestel-Clément-Sorine electromechanical model [Marchesseau 2012a, Marchesseau 2012b], first introduced by [Bestel 2001], was used (see Appendix B.7 for a description of the model and its parameters). Using SOFA allows efficient computation, modularity and interactivity which are strong assets for further therapy planning. This model includes available). We adjusted the initial set θ 0 after a trial and error approach on a healthy volunteer and then used the same initialization θ 0 for all studied cases (healthy and heart failure). Using the same θ 0 and range of parameter values for all cases makes the calibration automatic. In this respect, we have already presented consistent results for healthy controls as well as pathological cases in [START_REF] Marchesseau | [END_REF]].

B.3.2.2 ROUKF Personalization Algorithm

Once initialized, the parameters can be locally (regionally) estimated using a more advanced algorithm that deals with much larger dimensional systems. We use in this paper a type of Kalman filter which belongs to the sequential data assimilation methods.

Data Assimilation Principles

Data assimilation methods intend to estimate state variables of the considered system using measurements. State variables usually consist of the trajectory (positions, velocities ,...) but can in fact include model parameters or other dynamical variables such as regional strains for instance. A model is therefore represented by a dynamical system governing the state variable X that includes the parameters of the model:

Ẋ = A(X, t) X(0) = X 0 + ξ X (B.4)
where A is a non-linear operator, X 0 the a priori value for initial conditions and ξ X the uncertainties in the initial conditions. The available measurements are described through a linear observation operator H.

The observations are therefore modeled by the application of this operator:

Z = HX + χ (B.5)
where χ is the observation noise. Two types of assimilation methods exist: the sequential and the variational approaches. In the variational approach, a cost function is iteratively minimized after gathering measurements over a period of time. In the sequential approach, a new estimation of the state is performed each time a new measurement is available through statistical analysis. Both approaches have shown promising the parametric space of dimension p. Its rank is therefore p, much smaller that the dimension of the space d. We suppose that we can manipulate covariance matrices in the factorized form:

P = LU -1 L T (B.9)
where U is an invertible matrix of rank p and represents the main uncertainties of the system. A new derivation of the algorithm can therefore be performed without computing the full covariance matrix P, as shown in [Moireau 2011, Moireau 2008].

B.3.2.3 Observation Extraction from the Images

We suppose in this paper that the deformed meshes registered from the cine-MRI as explained in Sec.B.3.1.2 are the only available measurements. Several types of observations can then be extracted from these meshes. Since the registration method focuses only on the visible contours in the image, we cannot expect the trajectory of all the points in the myocardium to be accurate. We will therefore only consider the points on the endocardial surfaces. Three types of observations can be derived directly from this:

Volume curve First, for the calibration, the volume of the Left Ventricle is computed over time, giving the three required observations (V min , dV /dt min and dV /dt max ). These observations are easily computed using the topology of the endocardium as explained in [Marchesseau 2012a]. However it is not sufficient to recover the motion of the whole myocardium.

Positions on the Surface

The most complete observation include the positions of all the points on the endocardial surface. The estimation tries therefore to match the motion point-wise.

The advantage of what we call the Surface-Points Estimation is that the observation operator H is straightforward, filled with 1 values on each coordinate of the surface nodes and 0 values otherwise. The main drawback is that there is no guarantee that the registration tracks displacements of material points, but rather apparent displacements of the contours, inducing therefore errors on the observations (well known aperture problem in Computer Vision where several deformations would lead to the same visible contours). There exist observation operators [Chabiniok 2011] that filter this aperture issue taking into account distances between the surface contours. These methods give promising results but are not directly applicable for clinical context due to their computation times.

Regional Volumes

Alternatively, we suggest to use regional volumes as observations. To this end, we project the 17 AHA volumetric segments (Fig. B.6) on the LV endocardium surface, and calculate the volume formed by this surface and the barycenter of the LV endocardium over time. Not only does this method allow to capture the motion using Reduced Order UKF from Regional Volumes of the surface locally without smoothing tracking errors, but it also leads to a small number of observations to match, decreasing this way the computation time. Note that a possible drawback could be a loss of observability (or identifiability here) due to the reduction of the observation space. However, regional volumes are used by clinicians as a pathological index [Baxley 1971] and therefore may contain sufficient information on the investigated pathology. (Bottom) Example of Regional Volume.

B.3.2.4 Error Measurements

In order to evaluate the errors between the estimated motion from the registration and the simulated motion, several error measurements are used. First, the result of the calibration algorithm compares the simulated volume curve V simu (t) to the registered volume curve V obs (t) using the standard root mean square difference:

εV = Nsteps t=0 (V simu (t) -V obs (t)) 2 N steps (B.10)
where N steps is the number of images in the MRI sequence. Then, after personalizing the model from regional volumes, we compare the simulated regional volumes V simu i (t) to the registered ones V obs i (t) normalized by the area S i of the AHA surface (for an estimation of the distance between the surfaces) to obtain a mean square difference for each regional volumes that we average:

εRegV ol = 1 N N i=0 1 S i Nsteps t=0 V simu i (t) -V obs i (t) 2 N steps . (B.11)
It is also interesting to record the distance error between the regional volumes at end-systole where it is maximal:

E RegV ol = 1 N N i=0 V simu i (T max ) -V obs i (T max ) S i . (B.12) B.4. Results

131

These errors give a good indication on how well the simulated motion match the registered one on the left endocardium. For a more global index that takes all the surfaces into account (including epicardium and right endocardium), we define a symmetric distance on the contours Ω as:

ε C = 1 2 N points k=0 dist(Ω simu -P obs k ) N points + 1 2 N points k=0 dist(Ω obs -P simu k ) N points (B.13)
with a mean value over time noted εC and a value at end-diastole E C .

B.4 Results

B.4.1 Selection of Adequate Observations

In order to select the most appropriate type of observations (surface points vs LV endocardium regional volumes), we performed a synthetic study. We created a sequence of deformed meshes with a forward simulation of the model using different values of the contractility for each of the 6 zones of the left ventricle (combining the AHA segments as shown in Fig. B.7). We then estimated these 6 contractilities with the ROUKF personalization algorithm using either the surface points as observations or the regional volumes. The resulting deformed meshes after personalization are compared with the ground truth (simulated data) on Fig. B.8 where we can see that both methods give an excellent match. The resulting estimated parameters are given in Fig. B.9 for the two estimations and compared with the reference value. We can see that surprisingly the parameters are better estimated with the regional volumes. To understand why observing regional volumes leads to better parameters in this case than observing all the points of the surface, we evaluated the errors on these two types of observations (Surface Point error εC and Regional Volumes error εRegV ol ) while varying the contractility of one zone around its ground truth value. Fig. B.10 gives the resulting mean errors. We can see that the Regional Volumes error εRegV ol using Reduced Order UKF from Regional Volumes has only one minimum whereas the Surface Points error εC has a non-smooth variation leading to several local minima in which the estimated parameters may be "trapped" in the filtering process. Of course, observability is enhanced when using pointwise displacements, hence we conjecture that a finer tuning of the filtering method would overcome this difficulty, see [Imperiale 2011].

B.4.2 Detection of Infarcted Tissue

In this section, we tested the ability of a personalized model, using the regional volumes as observations, to detect regions of low contractility, for instance due to presence of scars following an infarct. In order to evaluate the detection power of the proposed personalization based on regional volumes, a set of tests were performed using synthetic data. Forward simulations with varying size of scar tissue (from 300 tetrahedra to 2000 tetrahedra out of the 65 000 tetrahedra for the full myocardium) were performed to create the observations. The relative contractility was set to 0.3 for the scar zone and 1 for the rest of the myocardium. Before using for the various sizes and noise amplitudes. We note that the noise does not influence the estimation of the parameter for scars zone larger than 1300 tetrahedra (as represented in Fig. B.12). Moreover, since a relative contractility of 1 represents the normal contractility for the myocardium, an estimated relative contractility too the mean error on the regional volumes is pretty low (0.7mm) which is smaller than a pixel size. However, the error on the total surface still needs to be improved including observations from the right ventricle and the epicardium.

B.4.4 Preliminary Specificity Study

Our database consisting of 8 healthy controls and 3 pathological cases selected for the Cardiac Resynchronization Therapy (with therefore a small ejection fraction and a long QRS) allows to draw preliminary conclusions and compare each pathological case to the controls. First, from the global volume curves, a calibration was performed to initialize 4 parameters. The resulting parameters for the healthy controls (box plots) and each pathological case are presented in Fig. B.19. We can see from this graph that the Computational Biophysical Model Personalized for case using Reduced Order UKF from Regional Volumes 1 (CBMP#1) and case 2 (CBMP#2) exhibit a similar behavior with a smaller global contractility than the healthy controls but normal bulk modulus, viscosity and peripheral resistance. On the other hand, the Computational Biophysical Model Personalized for the third case (CBMP#3) exhibits what seems to be a global contractility closer to the normal range, and a very high Bulk Modulus. Therefore, the model of the myocardium muscle for this case is stiffer than the model describing the healthy cases but the contractility is equivalent. The small ejection fraction can therefore be due for CBMP#3 to the high stiffness whereas it can be due for CBMP#1 and CBMP#2 to a small contractility. This difference could be explained by the fact that this patient has a different etiology than the first two cases (idiopathic cardiomyopathy) and, more importantly, has regurgitation at both the mitral and aortic valves. Furthermore, this patient presents significant atrial volumes, with an increased pulmonary artery pressure. This leads to a smaller inflow than outflow as shown in A comparison of the regional contractilities obtained from the personalization method that we propose leads to consistent observations (see Fig.B.21). Indeed, the personalized regional contractilities of the first two cases are much smaller on all the regions than the controls. Moreover, CBMP#3 exhibits regional contractilities closer to the contractilities of healthy hearts for most regions. Some of the estimated contractilities however are extreme (zone 7 has a very small contractility while 17 has a very high contractility) which could indicate local dysfunctions. The study of the covariance on the parameters showed that our confidence on the contractilities of zones 1 to 6 and zone 17 is 4 times smaller than on the other zones due to the boundary constraints on the base and the apex that we set. .21: Estimated regional contractilities for the 3 pathological cases compared to the healthy controls (box plots). using Reduced Order UKF from Regional Volumes

B.5 Discussion

B.5.1 Personalization from Regional Volumes

In this paper we proposed patient-specific modelling of the heart based on regional volumes. Tests performed on synthetic data proved that this choice of observations was relevant to estimate regional contractilities and detect scar tissue as long as the size of the scar zone is not too small. The advantage of such an approach is to smooth the input observations in order to match the motion of the left endocardium. This method proved to be robust and can be used at a large scale. Indeed, although the regional volumes were rather different for pathological cases and healthy volunteers, we managed to estimate the regional contractilities and have a good final match on the motion. Finally estimating the regional volumes gives precious indices on the synchronization of the AHA zones which could be further studied for therapy planning.

B.5.2 Analysis of Real Case Personalizations

Our database of 3 pathological cases and 8 volunteers allows to draw preliminary conclusions. First, the personalized computational biophysical models that intend to characterize these pathological cases exhibit different global contractility and/or stiffness. Our results on the calibration and on the regional personalization are consistent and show that CBMP#1 and CBMP#2 display a similar behavior (small estimated global and regional contractilities with normal estimated stiffness), while CBMP#3 is clearly different. First, CBMP#3 exhibits a reasonable global contractility parameter while showing a high stiffness which seems to explain its small ejection fraction and its small filling flow. Second, CMBP#3 estimated regional contractilities are for most of them similar to healthy cases with localized akinetic regions confirmed by the clinicians. In addition, this case, with a different etiology, was the one not responding to the CRT. Our personalization pipeline therefore gives results in agreement with the clinical interpretation of the images.

B.5.3 Limitations of the Personalization Pipeline

The proposed approach includes several processing and computational steps that are subject to uncertainties. First, we relied on image segmentation and tracking which are challenging due to the limited intensity contrast of the myocardium, and the aperture problem. Then, with the calibration, 4 global parameters could be assessed (contractility, viscosity, bulk modulus and peripheral resistance) and the remaining 10 (including other stiffness parameters, contraction and relaxation rates as well as other pressure related parameters) were set to their standard values which were chosen after a trial and error approach on one healthy case, or from the literature (see Appendix B.7 for the full list and values). This choice of parameters to calibrate was made from a complete sensitivity analysis [Marchesseau 2012a] which led to conclusive results on healthy and pathological cases [START_REF] Marchesseau | [END_REF]]. However, data on ventricular pressure would help the parameter estimation compared to using literature values that can be far from the actual patient condition. Second, boundary conditions play a significant role in the myocardium motion and hemodynamics, therefore they need to be assessed carefully. Our model includes a pericardium that constrains the outward radial motion and some elastic constraint that can be applied on the apex (which is either outside the field of view of the images, or barely moving) and on the basal area around the valves. However, the pre-load and the after-load were not personalized in this article which probably impacts the estimation of the parameters. For instance, modelling the regurgitation observed on the third pathological case could help in obtaining parameters closer to the actual intrinsic tissue properties of this patient. Finally, a time registration between the real data and the simulation is necessary and influences as well the results of the personalization.

B.5.4 Perspectives

Our method only intends to recover the motion of the left endocardium. To improve the results of the personalization, or estimate more parameters such as the regional stiffness for instance, one may include the regional volumes calculated from the epicardium and intend to recover the thickening and shrinking of the muscle or estimate the parameters related to the boundary conditions. Using measurements on the regurgitation would also improve the personalization since our heamodynamic model that represents the pressure constraint [Marchesseau 2012a] can include regurgitation by estimating the iso factor K iso . Alternatively, we could couple this heamodynamic constraint with a lumped-flow model as used by [Koon 2010]. Also, since the estimated parameters indicate a highly incompressible behavior, we will address this incompressible behavior in the future by adopting specific finite element procedures. Moreover, we could divide the endocardium in different zones than the standard AHA zones in order to obtain more precise results (for instance taking into account known scar zone, or dividing the actual AHA zones by 2). Applying this method on a larger database could also lead to better assessment of the parameters specificity of the personalized models. Finally, the proposed method could also help in testing the acute response to pacing by simulating it on the personalized model.

B.6 Conclusion

This paper proposed a personalization strategy made of a global calibration and the application of the Reduced Order Unscented Kalman Filter to estimate regional contractilities based on regional volumes using a complete model of the heart. The personalization of 3 pathological cases and 8 volunteers shows the robustness of this strategy and opens up possibilities to study the specificity of the estimated physiological parameters to the etiology and the response to CRT. The use of regional volumes allows to smooth registration errors and to find a good compromise between the physiological behavior of the model and the accuracy of the personalization.

B.7

The Bestel-Clément-Sorine electromechanical model

For the sake of completeness, we provide a brief description of the Bestel-Clément-Sorine electromechanical model [Bestel 2001] further improved by [Chapelle 2012b] that we used throughout this paper. This description has already been published in [Marchesseau 2012b].

The model is composed of a passive isotropic visco-hyperelastic component that accounts for the elasticity and the friction in the cardiac extracellular matrix surrounding the fibres, described as a MooneyRivlin material. The strain energy for a MooneyRivlin material is given as:

W e = c 1 ( Ī1 -3) + c 2 ( Ī2 -3) + K 2 (J -1) 2
where c 1 , c 2 are material parameters and K is the Bulk modulus. The quantities Ī1 and Ī2 are the isochoric invariants of the Cauchy-deformation tensor C.

In parallel, the stress along the cardiac fibre is composed of an active part (contraction in the sarcomere) and a passive part corresponding to the elastic bound (titin) between sarcomeres and Z-discs, having stress σ s = E s e s . The contractile component having stress tensor σ c , driven by the control variable u, has a viscous part to account for the energy dissipated in the sarcomere due to friction. This gives At the nanoscopic scale, the binding and unbinding process of the actin and myosin filaments in the sarcomere is described by Huxley's filament model [Huxley 1957]. Statistical mechanics allows to describe its behavior at the B.7. The Bestel-Clément-Sorine electromechanical model 143 macroscopic scale, resulting in a differential equation that controls the active stress τ c and the sarcomere stiffness k c :

σ c = τ c + µ ėc .
kc = -(| u | +α | ėc |)k c + n 0 k 0 | u | + τc = -(| u | +α | ėc |)τ c + ėc k c + n 0 σ 0 | u | +
where α is a constant related to the cross-bridge release due to a high contraction rate, k 0 and σ 0 are respectively the maximum stiffness and contraction. n 0 is a reduction factor that allows to take into account the Starling effect by which the maximum contraction depends on the fibre strain e c . The control variable u is derived from the electrical activation model and is a function of the free calcium concentration only. It is modeled using electrophysiological inputs such as depolarization times (T d ) and action potential durations (AP D) and depends on two parameters: k AT P the rate of the myosin ATPase activity controlling the contraction rate and k RS the rate of sarcoplasmic reticulum calcium re-uptake controlling the relaxation rate.

The ventricles are filled with blood coming from the atria and ejected through the arteries. A valve model gives relationships between the outward flow and the various pressures (ventricular, arteria and atria). The arteria pressure is modeled using a four-element Windkessel model [Stergiopulos 1999], that depends on four parameters: the peripheral resistance R p , the characteristic time τ , the characteristic resistance Z c and the total arteria inertance L.

The values of the mentioned parameters are given Table .B.1.

C.2 Introduction

Instead of the commonly used body planes (coronal, axial and sagital), the cardiac MR images are usually acquired along several oblique directions aligned with the structures of the heart. Imaging in these standard cardiac planes ensures efficient coverage of relevant cardiac territories and enables comparisons across modalities, thus enhancing patient care and cardiovascular research. Optimal cardiac planes depend on global positioning of the heart in the thorax. This is more vertical in young individuals and more diaphragmatic in elderly. Automatic recognition of this metadata is essential to appropriately select image processing algorithms, to group related slices into volumetric image stacks, to enable filtering of cases for a clinical study based on presence of particular views, to help with interpretation and visualisation by showing the most relevant acquisition planes, and in content based image retrieval for automatic description generation. Although this orientation information is sometimes encoded within two DICOM image tags: Series Description (0008,103E) and Protocol Name (0018,1030), it is not standardised, operator errors are frequently present, or this information is completely missing. Searching through large databases to manually cherrypick relevant views from the collections is therefore tedious. The main challenge for an image content-based automated cardiac plane recognition method is the variability of the thoracic cavity appearance. Different parts of organs can be visible even across the same acquisition planes between different patients. An excellent introduction to standard MR cardiac acquisition planes can be found in [START_REF] Taylor | [END_REF]]. These planes are often categorized into two groups -the short and the long axis planes. In this paper we learn to predict the five most commonly used cardiac planes acquired with steady-state free precession (SSFP) acquisition sequences. These are the short axis, 2-, 3-and 4-chamber and left ventricular outflow tract views. These five labels are the target of our learning algorithm. See C.1 for a visual overview.

C.2.1 Cardiac acquisition planes

The left ventricular short axis slices (C.1a) are parallel to the mitral valve ring. These are acquired regularly spaced from the cardiac base to the apex of the heart, often as a cine 3D+t stack. These views are excellent for reproducible volumetric measurements or radial cardiac motion analysis but their use is limited in atrio-ventricular interplay or valvular disease study.

The long axis acquisition planes include the 2-chamber, 3-chamber, and 4-chamber views (C.1b,C.1c,C.1d). These planes are used to visualize different regions of the left atrium, mitral valve apparatus, and left ventricle. The 3-chamber and left ventricular outflow tract (also known as the coronal oblique view) (C.1e) views provide visualization of the aortic root from two orthogonal planes. The 4-chamber view enables visualization of the tricuspid valve and right atrium.

C.2.2 Previous work

The previous work on cardiac view recognition has been concentrated mainly on real-time recognition of cardiac planes for echography [Otey 2006, Park 2007a, Park 2007b[START_REF] Beymer | [END_REF]. There exists some work on magnetic resonance [Zhou 2012, Margeta 2014a, Shaker 2014]. The common methods are based on dynamic active shape models [START_REF] Beymer | [END_REF]], require to train part detectors [Park 2007b] or landmark detectors [Zhou 2012]. Therefore any new view will require these extra annotations to be made. [Otey 2006] avoid this limitation by training an ultrasound cardiac view classifier using gradient based image features. [Margeta 2014a] trained classification forests to predict the cardiac planes directly from cardiac MR image miniatures. [Shaker 2014] recently proposed a technique based on autoencoders. They learn image representations in an unsupervised fashion and use this representation to distinguish between two cardiac views. We will later show that our proposed method reaches a very favourable performance and surpasses the previously introduced [Margeta 2014a]. We report the results on an open cardiac dataset which simplifies direct comparison of these different techniques in the future.

The state of the art in image recognition has been heavily influenced by the seminal works of [Krizhevsky 2012] and [Ciresan 2012] using convolutional neural networks (CNNs). [Krizhevsky 2012] trained a large (60 million parameters) CNN on a massive dataset consisting of 1,2 million images and 1000 classes [Russakovsky 2014].

They employed two major improvements: Rectified linear unit nonlinearity to improve convergence, and Dropout [START_REF] Hinton | [END_REF]]. The Dropout means that during the training phase the output of each hidden neuron is dropped out (set to 0) with a certain probability p. The dropped out neurons do not contribute to the forward pass and are skipped in the back-propagation. For every training batch a different set of neurons is dropped out i.e. a different network architecture is sampled. At test time, however, all neurons are used and their output responses are multiplied by 1 -p to compensate for the architecture difference. [START_REF] Hinton | [END_REF] showed that this strategy helps to reduce complex coadaptations of the neurons (as neurons learn not to rely on single neuron responses from the preceding layer) and to reduce overfitting.

Training a large network from scratch without a large number of samples still remains a challenging problem. A trained CNN can be adapted to a new domain by reusing already trained hidden layers of the network, though. It has been shown e.g. by [Razavian 2014] that the classification layer of the neural net can be stripped, and the hidden layers can serve as excellent image descriptors for a variety of computer vision tasks (such as for photography style recognition by [Karayev 2013]). Alternatively, the prediction layer model can be replaced by a new one and the network parameters can be finetuned through backpropagation. In this paper we consider all of these approaches (training a network from scratch, reusing a hidden layer features from a network trained on another problem, and finetuning of a pretrained network) for using a CNN in cardiac view recognition and compare it with with the prior methods based on random forests using pixel intensities from image miniatures or image plane normal vectors as features.

In this paper we compare the three groups of methods for automatic cardiac acquisition plane recognition. The first one is based on DICOM meta information, and the other two completely ignore the DICOM tags and learn to recognize cardiac views directly from image intensities. In C.3.1 we first present the recognition method using DICOM-derived features (the image plane orientation vectors). Here, we train a random forest classifier using these 3-dimensional feature vectors. We recall (C.3.2) the previously presented random forest based method [Margeta 2014a]. The newly proposed third path using convolutional neural networks is described in C.3.3. The later two learn to recognize cardiac views from image content without using any DICOM meta information. In C.4 we compare these approaches. Finally, in C.5 we present and discuss our results.

C.2.3 Contributions

The main contribution of our paper is that a good cardiac view recognizer (reaching state of the art performance) can be efficiently trained end to end without designing features using a convolutional neural network. This is possible by finetuning parameters of a network previously trained for a large scale image recognition problem.

• We achieve state of the art performance in cardiac view recognition by learning an end to end system from raw image intensities using CNNs

• This is one of the first papers to demonstrate the value of features extracted from medical images using CNNs trained on a large scale visual object recognition dataset

• We show that finetuning weights of CNN pretrained on an object recognition dataset is a good strategy that helps to improve performance and speed up the network training

• We show that the CNNs can be applied to smaller datasets (a common problem in medical imaging) thanks to careful network initialisation and dataset augmentation

C.3 Methods

A ground truth target label (2CH, 3CH, 4CH, LVOT or SAX) was assigned to each image in our training set by an expert in cardiac imaging. We use it only in the training phase as a target to train and to validate the view recognizers. At the testing phase we will predict this label from the features. In the paper, we compare two types of methods. The first one is using 3-dimensional feature vectors (image plane normal vector) computed from the DICOM orientation tag (C.3.1) and a random forest classifier which is trained to predict cardiac views using this vector. Algorithms in DICOM tag removal after an incorrectly configured anonymisation procedure, when parsing images from clinical journals or using other image formats. In these cases we have to rely on recognition methods using exclusively the image content. In the next two sections, we present two such methods. One that is based on previous work using decision forests and image miniatures [Margeta 2014a] and the other one is a new approach using convolutional neural networks. We learn to predict the cardiac views from 2D image slices individually, rather than using 2D + t, 3D or 3D + t volumes. This decision makes our methods more flexible and applicable also to view recognition scenarios when only 2D images are present e.g. when parsing clinical journals or images from the web. This simple method [Margeta 2014a] was posed as a standard image recognition framework where features are extracted from the images and are fed to a classifier (C.3), in this case a classification forest. Classification forest [START_REF] Breiman | Random forests-random features[END_REF]] is an ensemble learning method that constructs a set of randomized binary decision trees. Its advantage is computational efficiency and automatic selection of relevant features. At training time the tree structure is optimized by recursively partitioning the data points into the left or the right branches such that points having the same label are grouped together and points with different labels are put apart. Each node of the tree sees only a random subset of the feature set and greedily picks the single most discriminative feature. This helps to make the trees in the forest different from each other and leads to better generalization. At test time features chosen at the training stage are extracted and the images are passed through the trees and reach a set of leaves. Class distributions of all reached leaves are averaged across the forest and the most probable label is selected as the image view. See [Criminisi 2011b] for a detailed discussion on decision forests.

C.3.2.1 Using image miniatures

The radiological images are mostly acquired with the object of interest in the image center and some rough alignment of structures can be expected (See C.4). Image intensity samples at fixed positions (without extra registration) therefore provide strong cues about the position of different tissues (e.g. dark lungs in the top left corner or bright cavity in the center).

The central square from each image is extracted, resampled to 128x128 pixels and linearly rescaled to range between 0 and 255. We subsample the cropped centers to two fixed sizes (20x20 and 40x40 pixels). In addition we divide the image into non-overlapping 4x4 tiles and compute intensity minimum and maximum for each of these tiles (C.5). This creates a set of pooled image miniatures (32x32 pixels each).

The pooling adds some invariance to small image translations and rotations (whose effect is within the tile size). The pixel values at random positions of these miniature channels are then used directly as features. At each node of the tree 64 random tile locations across all miniatures are tested and the best threshold on this value is selected to divide the data points into the left or the right partitions until not less than 2 points are left in each leaf. We train 160 such trees. This method is shown in the evaluation as Miniatures + forest.

C.3.2.2 Augmenting the dataset with jittering

While the object of interest is in general placed at the image center, differences between various imaging centers and positioning of the heart on the image exist. The proposed miniature features are not fully invariant to these changes. To account for this the training set was augmented with extra images created by transforming the originals. These were translations (all shifts in x and y between -10 and 10 pixels for a 5x5 grid), but also scale changes (1-1.4 zoom factor with 8 steps while keeping the same image size) and in-plane rotations around the image center (angles between -10 and 10 degrees with 20 steps). The augmented images were resampled with linear interpolation. Note that the extra expense of dataset augmentation is present mainly at the training time. The test time remained almost unaffected except that now a deeper forest could be learnt. The benefit of dataset augmentation is clear, yielding a solid 12.14% gain in the F1 score (F 1 = 2(precision.recall)/(precision + recall)). Results using this augmented dataset are presented in the evaluation as Augmented miniatures + forest. Compared to the preliminary study in [Margeta 2014a], we will test this method on a bigger dataset.

C.3.3 Convolutional neural networks for view recognition

To improve the performance of the forest-based method, we turn to the state of the art in image recognition -the convolutional neural networks. Their principle is rather simple -the input image is convolved with a bank of filters (conv) whose response are fed through a nonlinear function e.g. a Rectified Linear Unit (ReLU) Acquisition Plane Recognition f (x) = max(0, x). The responses are locally aggregated with max-pooling. The output of the max-pooling creates a new set of image channels which are then fed through another layer of convolutions and nonlinearities. Finally, the fully connected layers (fc) are stacked and connected to a multiway soft-max in order to predict the target label. The role of the max-pooling is similar to the forest based method i.e. to aggregate local responses and to allow some invariance of the input to small transformations. The parameters of the network such as weights of the convolutional filters are optimized through backpropagation. Using the stochastic gradient descent the average batch prediction loss (soft-max) is decreased.

We use the widely adopted neural network architecture (see C.6) described by [Krizhevsky 2012] as implemented in the Caffe framework [Jia 2014] under the bvlc_reference_caffenet acronym (CaffeNet in short). The CaffeNet implementation differs from Krizhevsky's AlexNet in the order of local response normalisation (LRN) and max-pooling operations. Not only AlexNet/CaffeNet is the winning entry of the ImageNet LSVRC 2012 competition, the weights and definitions of this network are publicly available thus reducing computational time needed for the training and improve reproducibility of this work. We initialized the network weights from a pretrained CaffeNet. We then replaced the last 1000-dimensional logistic regression layer (previous fc8) with a 5-dimensional for 5 cardiac views. Then we finetuned the network with our data. We also extracted features from the last 4096-dimensional fully connected layer fc7 (in dark gray) from both CaffeNet and our finetuned Car-dioViewNet and used them with a linear SVM and a decision forest classifier. We achieve the best performance with our finetuned network directly. ReLU -rectified linear unit, max -max pooling, LRN -local response normalisation, conv-convolutional layer, fc -fully connected layer.

C.3.3.1 Classifying cardiac views using CNN features optimised for visual recognition

Similarly to the work of [Karayev 2013] for photography style recognition we use CNN features from a fully connected layer of a network pretrained on the ImageNet LSVRC [Russakovsky 2014] dataset. The fully connected layer (in our case fc7see C.6) helps us to describe the cardiac images with 4096-dimensional features.

Before putting the cardiac images through the network, simple preprocessing is done. The cardiac images are resized to 256x256 squares regardless of their input dimensions. Since the CaffeNet takes RGB images as input we simply replicate the 2D cardiac slices into each colour channel. We compute a pixel-wise mean image for our cardiac view dataset and subtract it from all training and testing images prior to entering the CNN. This centers image intensities around zero and serves as a very simple intensity normalisation. As we later found in our case, the average image is almost uniform and a scalar value could be subtracted instead. The central (227x227x3) crop of this image is then fed forward through the network. We then use these CaffeNet fc7 features with a linear SVM classifier to predict the cardiac views. We ran crossvalidation on the training subset folds to maximise the prediction precision. This helped us to choose the penalty parameter C of the SVM classifier from standard set of values [0.1, 1, 10, 100, 1000] as C = 1. We report results of this method as CaffeNet fc7 + SVM. Similarly, we trained a classification forest (with 64 features tested per node and 1000 trees) using these features (instead of image miniatures) and report these results as CaffeNet fc7 + forest. These features were adapted to a general object recognition task and come from a CNN that never saw a cardiac MR image to optimise its weights. As we will show in C.1, this already performs quite well for the cardiac view recognition. In the following we will show how we can further improve performance by adapting the CNN weights to the medical imaging domain.

C.3.3.2 CardioViewNet architecture and finetuning the visual recognition CNN

In practice, many examples are often needed to train a large capacity neural network. However, by starting from the weights of a pretrained network we can just finetune the network parameters with new data and adapt it to the new target domain. Here we use the pretrained CaffeNet [Jia 2014] and replace the last 1000-class multinomial regression layer with a 5-class one (See C.6). The net is finetuned with stochastic gradient descent with higher learning rate (10 -2 ) for parameters in the newly added layer and smaller (10 -3 ) in the rest of the network. We use momentum of 0.9 in the stochastic gradient descent optimiser, and a small weight decay at each iteration (10 -4 ). The batch size used in each iteration was 32 and the step size is kept constant for the whole training. A set of resized 256x256x3 images is used for training. At each iteration a batch of 32 random (not necessarily central) 227x227x3 crops is extracted from the resized cardiac slices and is fed forward through the network.

Compared to the implementation of the forest based method where all augmented images were precomputed, the translations are cheaply generated on the fly at each iteration. Already after 3000 iterations the prediction error on the validation dataset reaches a plateau and further improvements are only marginal, see C.7). We stop the optimization at 8000 iterations and pick this model in our experiments as it yields the best performance. To reduce overfitting, we use the Dropout strategy [START_REF] Hinton | [END_REF]] in the fully connected layers fc6 and fc7 with probability of dropping target (similarly to the CardioViewNet). The second network architecture is the CardioViewNet (CardioViewNet from scratch). We found the choice of the learning rate (10 -3 for CardioViewNet and 10 -5 for LeNet-5, both using batch sizes of 32) and good random initialisation to be crucial to avoid divergence. We initialise the weights with the procedure described by [He 2015] that is well suited for networks with the ReLU nonlinearity and choose the learning rate by trial and error, i.e. reducing the learning rate until the network starts to reduce the loss without diverging. We ran a randomized 10-fold cross validation by taking a random subset of 90% of the patients (rather than image slices) for training and used remaining 10% for validation. The patient splits guarantee that repeated acquisitions from the same patient that are occasionally present in the dataset never appear in both the training and the validation set and do not bias our results. Classification accuracy is not a Table C.1: Evaluation of the algorithms in the two groups of algorithms, we highlight in bold the best performance from each group. References to relevant sections in the paper with more details on each algorithm in the parentheses. We computed average of individual view F1 scores, precisions and recalls for each fold (except for the underrepresented LVOT) for the two datasets. Here we display means and standard deviations of these average scores across all 10 folds. The prediction using classification forests on the DICOM orientation vector is the best performing method. However from purely image based methods, the finetuned convolutional network CardioViewNet outperforms the rest. good measure for imbalanced datasets as performance on the dominant class (i.e. short axis) can obfuscate the true performance. Therefore in this paper we report means and standard deviations of average (averaging is done across the classes) precisions, recalls and F1 scores. In the context of content based retrieval, these measures can be interpreted as following: The precision (also known as positive predictive value or false positive rate is defined as TP / (TP + FP )) is the fraction of relevant images (having the same view as the query view) out of all returned images. The recall (also known as sensitivity or true positive rate is defined as TP / (TP + FN ) ) is the probability of retrieving a relevant image out of all existing relevant images. Where TP is the number of true positives, FP the number of true negatives, and FN the number of false negatives. The F1 score is the harmonic mean of the precision and the recall.

C.4 Validation

To study the robustness of the presented algorithms against the dataset bias, we trained recognizer models from different folds trained exclusively on the DE-TERMINE dataset ( [Kadish 2009]) and tested them on a completely independent dataset -the STACOM motion tracking challenge [Tobon-Gomez 2013] (KCL in short) containing slices from 15 patients (2CH:15, 4CH:15, SAX:207). The KCL dataset consists of healthy volunteers and the images are in general of higher and more uniform quality and with more consistently chosen regions of interest. This allows us to evaluate performance on the 3 cardiac views present. We invite the interested readers to look at this open access dataset through the Cardiac Atlas Project website [Fonseca 2011].

C.5. Results and discussion

159

C.5 Results and discussion

Here we present results for the method using DICOM based prediction and methods using image content. The mean average F1 scores, precisions and recalls are summarized in C.1 and total confusion matrices for the two best methods from each family (DICOM based and Image based) can be seen in C.9. We confirm findings from the previous work that cardiac views can be predicted from image plane normals and serve as a prior [Zhou 2012] but require presence of the relevant DICOM tag. The larger DETERMINE dataset turned out to be more challenging for the miniature based method and it performed significantly worse than in the results published previously. On the contrary the performance of the CaffeNet features for cardiac view description is quite remarkable. These were not trained for cardiac MR images yet they perform better than most methods with handcrafted features. They most likely encode local texture statistics which helps with the prediction. Adding some texture channels to the miniature method could therefore improve the performance. The quality of predictions using the finetuned CardioViewNet is almost on par with the approach using DICOM derived information and significantly outperforms the previous forest based method using image miniatures [Margeta 2014a] while not requiring to train any extra landmark detectors as in [Zhou 2012]. As features extracted from the CardioViewNet do a good job even when used with external classifiers, they could be used to learn extra view recognizers without additional finetuning on these. In C.10 we present examples of some of the least confident correct predictions using the finetuned CardioViewNet. It is important to note that the softmax output of the neural network not only returns the label but also some measure of confidence in the prediction. Similarly, the incorrectly classified images (C.11) often belong to views more difficult to recognise for a human observer and the second best prediction is the correct result. When training from scratch, the performance does not seem to improve beyond 10000 iterations and we pause the backpropagation there. Although the performance is much lower than for the finetuned networks, the network learns to predict the cardiac views. For the LeNet-5 model trained from scratch, the performance closely follows the CaffeNet fc7 + SVM model. We did not observe any benefit of using forests (at least when using orthogonal splits) over the linear SVM when using features from the convolutional nets and the forest performs in general worse. The performance of the CardioViewNet trained from scratch is better than using the fully connected layer features (fc7 ) from CaffeNet but the training takes significantly longer to obtain.

The predictions on the KCL dataset are naturally slightly worse as some differences between the studies still exist. We have observed that test time oversampling (averaging predictions of the central and corner patches and their horizontally flipped versions) improves the scores for this dataset although it does not improve the DETERMINE dataset predictions. This might indicate that a better thought oversampling, image normalization or dataset augmentation strategies might further improve the cross-dataset generalisation.

C.6 Conclusion

Convolutional neural networks and features extracted from them seem to work remarkably well for medical images. As large datasets to train complex models are often not available, retargeting the domain of a previously trained model by finetuning can help. This speeds up the learning process and achieves better predictions. Even models trained for general object recognition might be a great baseline to start. In our case doing network surgery and finetuning the pretrained model allowed us to make significant progress in cardiac view recognition from image content without Similarly for the 4CH view, the right atrium is missing and one can already see parts of the outflow tract branching out of the left ventricle typical for a 3CH view. The 3CH view is captured with a detailed region of interest not very common in the dataset. Extra data augmentation could probably help to fix this case. Finally, the LVOT views are severely under-represented and can be confused with basal short axis views. Note that for all of these cases, the correct prediction is in the top 2 most likely views.

[Cheng 1999] Jianhua Cheng, Kaichiro Kamiya, Weiran Liu, Yukiomi Tsuji, Junji Toyama and Itsuo Kodama. Heterogeneous distribution of the two components of delayed rectifier K+ current: a potential mechanism of the proarrhythmic effects of methanesulfonanilideclass III agents. Cardiovascular Research, vol. 43, no. 1, pages 135-147, 1999

Figure 1

 1 Figure 1.2: Manuscript Organization.

Figure 2

 2 Figure 2.1: a) Heart anatomy and chambers. b) Layers of the heart wall. From en.wikipedia.org/wiki/Heart.

Figure 2

 2 Figure 2.2: Human ventricular fibers from in vivo Diffusion Tensor Imaging (DTI). Fibers are color-coded with the local helix angle α. Courtesy of team.inria.fr/ asclepios.

Figure 2

 2 Figure 2.3: Heart's conduction system. From en.wikipedia.org/wiki/Heart.

Figure 2 . 4 :

 24 Figure 2.4: (Left) Simplified schematic of ionic flow through the cell membrane. (Center) Evolution of ionic conductances. (Right) Action potential (mV). Courtesy of [Talbot 2014]

Figure 2

 2 Figure 2.5: Three-dimensional reconstruction of infarct border zone. (Left) Shortaxis slice cut 4mm of the left ventricular apex. Collagen is indicated in white/yellow and myocytes by red/brown. (Right) Reconstructed infarct border zone. The dimensions of the image volume are 2.99×2.68×0.70mm 3 . Courtesy of [Rutherford 2012]

  ]: -Non-fluoroscopic localization of the ablation catheter within the heart -Display of electrogram characteristics (i.e. activation time or voltage) w.r.t. anatomic position -Integration of electroanatomic information with 3D images of the heart obtained from computed tomography (CT), magnetic resonance imaging (MRI) or other imaging technique.

Figure 2 . 6 :

 26 Figure 2.6: Example of matching histological data in chronic infarcts against DE-MDCT (A and B) and DE-MRI (D and E). The peri-infarct zone (PIZ) is visualized between tow black arrows by intermediate signal intensity (white circle) in DE-MDCT images. In D and E, DE infarct scar (white arrows) and the PIZ (white circle) show different signal intensity in DE-MRI. (C and F) Masson trichrome stain depicts viable myocardium in red ( * ) from non-viable tissue in blue. At higher magnification (F) the islands of viable myocytes (red) within the scar tissue are visualized, showing the heterogeneity of the PIZ. Image from [Schuleri 2009].

Figure 2

 2 Figure 2.7: CARTO navigation system with (Left) fluoroscopic images, (Center) CARTO contact mapping catheters and (Right) the electroanatomical mapping. (Images from CHU, Bordeaux)

Figure 2

 2 Figure 2.8: PentaRay multi-spline recording catheter. (Image from BioSense Webster Inc.)

Figure 2 . 9 :

 29 Figure 2.9: EnSite navigation system with (Left) fluoroscopic images, (Center) the deflated and inflated balloon and (Right) the electroanatomical mapping. (Images from KCL, London)

  Figure 2.10.

Figure 2 .

 2 Figure 2.10: Sample unipolar (top and middle) and bipolar (bottom) intracardiac electrograms obtained from the CARTO navigation system.

  zu 2 (1-u) τ in represents the inward currents, J in , which raise the action potential voltage (primarily N a + and Ca 2 +) u τout represents the outward currents, J out , that decrease the action potential voltage (mainly K + ), describing repolarization J stim is the stimulation current at the pacing location τ in , τ out , τ open , τ close have units of seconds D is a diffusion tensor that controls the diffusion term in the model This model incorporates both action potential duration (APD) and conduction velocity (CV) restitution effects, and the restitution curves can be written in an analytical formulation. APD restitution. It is an electrophysiological property of the cardiac tissue and defines the adaptation of APD as a function of the heart rate. Its slope has a heterogeneous spatial distribution. The APD restitution curve (APD-RC) defines the relationship between the diastolic interval (DI) of one cardiac cycle and the APD of the subsequent cardiac cycle. The slope of these RCs is controlled by a model parameter τ open of the MS model and depicts the APD heterogeneity present at multiple heart rates. APD-RC for MS model is explicitly derived as:

Figure 2 .

 2 Figure 2.11: Classification: Training data and tree training. a) Input data points. The ground truth label of training points is denoted with different colors. Grey circles indicate unlabeled, previously unseen test data. b) A binary classification tree. During training a set of labeled training points {v} is used to optimize the parameters of the tree. In a classification tree the entropy of the class distributions associated with different nodes decreases (the confidence increases) when going from the root towards the leaves. Courtesy of [Criminisi 2011a]

Figure 2 . 12 :

 212 Figure 2.12: Effect of forest size in classification. Increasing the forest size T produces smoother class posteriors. Courtesy of [Criminisi 2011a]

Figure 2 . 13 :

 213 Figure 2.13: Classification forest testing. During testing the same unlabeled input data v is pushed through each component tree. At each internal node a test is applied and the data point sent to the appropriate child. The process is repeated until a leaf is reached. At the leaf the stored posterior p t (c|v) is read off. The forest class posterior p(c|v) is simply the average of all tree posteriors. Courtesy of [Criminisi 2011a]

Figure 3 . 1 :

 31 Figure 3.1: Personalized computer modelling process. Upper panel: (A) highresolution contrast-enhanced MRI scar images; (B) whole heart model segmented from 3D steady-state free precession (SSFP) MRI with scar (core and gray zone) in violet; (C) low voltage areas from electroanatomical mapping. Lower panel: (D) model personalization and in silico VT stimulation study procedure workflow.

Figure 3

 3 Figure 3.2: Box-whiskers plots of spatial heterogeneity of apparent conductivity (A); maximum APD-RC slope (B); and maximum APD (C) for each of the patients. Box = median and quartile ranges; whiskers = total range. Data range values for each of the boxplots are shown at the bottom of each box-plot. Colour coding shows the range varies from the largest (darker brown) and the smallest (lighter brown) for each parameter.

Figure 3

 3 Figure 3.3: LV polar plot representation of the spatial distribution of the maximum slope for APD-RCs, maximum asymptotic APD and AC. Scars identified by DE-MRI with its border zone (highlighted by gray contours) and clinically observed exit points (highlighted by white contours, fang-shaped due to unfolding of 3D volume surface to 2D polar plot) during re-entrant VT are overlaid to the polar plots.

Figure 3 . 4 :

 34 Figure 3.4: Percentages of LV spatial domain with high gradient (≥ 90 percentile differences between neighbouring regions) in terms of AC, maximum asymptotic APD and APD-RC slope for each patient.

Figure 3 . 6 :

 36 Figure 3.6: VT isochrones of induced re-entrant VT during clinical VT stimulation study (left panel) and during in silico VT stimulation study (right panel). The arrows point towards the directions of propagation. Scars (core and gray zone) are superimposed to the LV shell shown in darkened regions.

Figure 3

 3 Figure 3.7: Correlation of predicted exits points with structural heterogeneity and functional heterogeneity. A = Patient 1. B = Patient 2. The left panel shows the isochrones for the three predicted VT circuits with fangs of white lines denoting the exit points. The middle panel demonstrates the LV scar distribution with black region denoting healthy myocardium, white region denoting the scar core and the gray region denoting the gray zone/scar boarder zone with overlying composite exit points (fanged green lines). The right panel demonstrates the electrical properties with overlying composite exit points (fanged green lines).

  Based on: [Cabrera-Lozoya 2014] Cabrera-Lozoya, Rocio, et al. "Confidence-based Training for Clinical Data Uncertainty in Image-based Prediction of Cardiac Ablation Targets." bigMCV Workshop MICCAI 2014. 2014. [Cabrera-Lozoya 2015d] Cabrera-Lozoya, Rocio, et al. "Local late gadolinium enhancement features to identify the electrophysiological substrate of post-infarction ventricular tachycardia: a machine learning approach." Journal of Cardiovascular Magnetic Resonance 1.17 (2015): 1-2.

  EP Mapping. The CARTO mapping system (Biosense Webster) enables the 3D localization of the catheter tip and provides the distribution of electrophysiological signals on cardiac surfaces. Contact mapping was achieved in sinus rhythm on the endocardium (trans-septal approach) and the epicardium (sub-xiphoid approach) with a multi-spline catheter (PentaRay, Biosense Webster), screenshots of the resulting study are shown in Figure4.1.

Figure 4 .

 4 Figure 4.1: Screenshots of electroanatomical mapping data for Patients 1 (left), 2 (center) and 3 (right). The mesh's colormap indicates local activation times (LAT) values, crosses and spheres represent electroanatomical points.

Figure 4 .

 4 Figure 4.2: Electrical readings from two electroanatomical points depicting a healthy signal (blue) and a LAVA (red). Signals correspond to: (top) the bipolar signal obtained from the two distal most electrodes in the catheter, (middle) the unipolar potential from the distal most electrode and (bottom) the bipolar signal obtained from the difference of the two most proximal electrodes in the catheter.

Figure 4 .

 4 Figure 4.3: Spherical representation of the catheter's hearing range (blue) for a given location. Scar regions are shown in black.

  Figure 4.6.

Figure 4 . 4 :

 44 Figure 4.4: [left] Anatomical model and electro-anatomical points (blue: healthy, red: LAVA) in scar regions (black). [right] CARTO reconstruction of the endocardial cavity, with activation times.

Figure 4 .

 4 Figure 4.5: Location of two EAPs and temporal position variation described by a red ellipsoid (above: low displacement, below: high displacement).

Figure 4 . 6 :

 46 Figure 4.6: Posterior endocardial maps for uncertainty attributes [mm] (Left: projection distance, Right: major ellipsoid axes).

Figure 4 .Figure 4

 44 Figure 4.7: Precision-recall and ROC curves for patients 1-3 after nested crossvalidation with precision optimisation. Line colors represent each fold (curves with AUC = 1 are aligned with axes) and the dotted line represents the average curve for all folds.

. 4 .Figure 4

 44 Figure 4.8: Precision-recall and ROC curves for patients 1-3 after nested crossvalidation with precision optimisation and temporal displacement weighting. Line colors represent each fold (curves with AUC = 1 are aligned with axes) and the dotted line represents the average curve for all folds.

Figure 4 . 9 :

 49 Figure 4.9: Bar plots for classification performance scores using sample uncertainty weights

Figure 4 .

 4 Figure 4.10: [left] LAVA regions from ground truth. [center] LAVA regions from predictions. [right] Prediction errors. (Color coding by classification confidence)

Figure 5 . 1 :

 51 Figure 5.1: Heart mesh depicting healthy myocardium (yellow), scar (black) and greyzone (grey) tissues derived from DE-MRI segmentation.

Figure 5

 5 Figure 5.2: QRS window detection in six surface ECGs (left column: limb leads I, II, III. right column: augmented limb leads: aVL, aVR, aVF). Raw signals (blue), filtered signals (red), candidate QRS windows (green) and final QRS windows obtained through majority voting (black). QRS windows have been given an offset for better visualisation.

  Unipolar electrograms are obtained by positioning the exploring electrode in the heart and the second electrode (the indifferent electrode) away from the heart such that it has little or no cardiac signal[Stevenson 2005].

Figure 5

 5 Figure 5.5: (Left) Points depicting the 1472 electrode locations used for electrogram recording during the EP study of Patient 1. (Right) Electrode locations where the potentials were computed. M1 (blue), M2 (green), M3 (yellow), M4 (red).

Figure 5 .

 5 Figure 5.6 provides further insight in the distribution of the clinical annotations of the electrode locations. It can also be seen that LAVA annotated signals have a higher tendency to remain close to the border zone, nonetheless, healthy signals can also be found in these regions.

Figure 5 . 6 :

 56 Figure 5.6: (Left) Image-driven personalized model and location of LAVA (red) and non-LAVA (blue) annotated signals on Patient 1. (Right) Boxplots representing the distance of the annotated signals to the border zone. The same relation of distance to border zone and electrogram annotation were found for the other four patients.

Figure 5 . 7 :

 57 Figure 5.7: Simulated and clinical non-LAVA EGM with electrode locations in the healthy myocardium

Figures 5

 5 Figures 5.7, 5.8, 5.9 and 5.10 show sample EGM simulation signals (unipolar M1 and M2 signals and the corresponding M1-M2 bipolar signal), the location of the electrodes with respect to the cardiac geometry and the clinical CARTO signals read at these locations. Blue bipolar signals correspond to those with non-LAVA clinical tags, whereas the red ones correspond to LAVA.

Figure 5 . 8 :

 58 Figure 5.8: Simulated and clinical non-LAVA EGM with electrode locations in the healthy myocardium

Figure 5 . 9 :

 59 Figure 5.9: Simulated and clinical LAVA EGM with electrode locations in the scar and grey zone

Figure 5 .

 5 Figure 5.10: Simulated and clinical LAVA EGM with electrode locations in the scar and grey zone

Figure 5 .

 5 Figure 5.12: Boxplots showing the distribution of the number of inflection points of the LAVA (red) and non-LAVA (blue) class in simulated EGMs for the five patients. KS tests demonstrated significant difference between distributions.

Figure 5 .

 5 Figure 5.13: Boxplots showing the distribution of the mean slope of the LAVA (red) and non-LAVA (blue) class in simulated EGMs for the five patients. KS tests demonstrated significant difference between distributions.

Figure 6 .

 6 Figure 6.1: (Left) Anterior and (Right) posterior resulting regions for classification result display after endocardial surface division.

Figure 6 . 2 :

 62 Figure 6.2: Image-based Learning. (Left) LAVA regions from ground truth. (Center) Predicted LAVA regions with confidence level color-coding. (Right) Prediction errors with confidence level color-coding.

Figure 6 . 3 :

 63 Figure 6.3: Simulation-based Learning. (Left) LAVA regions from ground truth. (Center) Predicted LAVA regions with confidence level color-coding. (Right) Prediction errors with confidence level color-coding.

Figure 6 .

 6 Figure 6.4 show the merged results for image and simulation features by selecting the predictions with the highest confidence for Patient 1. The maximum confidence fusion scheme had overall high scores across patients. A slight increase in the classification accuracy (95.0%), LAVA-specificity (99.2%) and LAVA-precision (95.4%) is obtained with respect to the results yielded by the simulation-based classification. This happens at the expense of LAVA-sensitivity (63.0%). Nevertheless, the value is still considerably higher than that obtained by image-based classification alone.

Figure 6 . 4 :

 64 Figure 6.4: Maximum Confidence Fusion. (Left) LAVA regions from ground truth. (Center) Predicted LAVA regions with confidence level color-coding. (Right) Prediction errors with confidence level color-coding.

Figure 6 .

 6 Figure 6.5: Image-based, Simulation-Enhanced. (Left) LAVA regions from ground truth. (Center) Predicted LAVA regions with confidence level color-coding. (Right) Prediction errors with confidence level color-coding.

Figure 6 . 6 :

 66 Figure 6.6: Feature Set Fusion. (Left) LAVA regions from ground truth. (Center) Predicted LAVA regions with confidence level color-coding. (Right) Prediction errors with confidence level color-coding.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 7.2 Methodological and Clinical Perspectives . . . . . . . . . . . 101 7.2.1 Methodological Perspectives . . . . . . . . . . . . . . . . . . . 101 7.2.2 Clinical Perspectives . . . . . . . . . . . . . . . . . . . . . . . 102

Figure 7

 7 Figure 7.1: Python-based graphical user interface.

  As shown in representative MAPs of LV base of individual transgenic LQT2 and LMC rabbits (Figure A.1B) and in averaged bull's eye plots of AP D 75 in all LV regions (Figure A.1C), LQT2 and E4031-treated rabbits had significantly longer APDs than LMC. Moreover, spatial dispersion of APD was significantly greater in LQT2 and E4031 than in LMC (Figure A.1D). A similar difference observed in surface ECG with significantly longer QT intervals in LQT2 and E4031 than in LMC (Figure A.1E).

  Figure A.2A). Global peak systolic longitudinal velocities (Vz) did not differ. Regional systolic Vz, however, were reduced only in LQT2 in anteroseptal segments (Figure A.2B). Ejection fractions were normal in all groups (Figure A.2C). LV end diastolic and end systolic volumina (LVEDV, LVESV) were lower in the smaller LQT2 rabbits (Figure A.2D). After body weight-correction, however, no differences were found between groups (Figure A.2D).

  .), primarily located in mid (Figure A.5C).

Figure A. 1 :

 1 Figure A.1: Spatial Differences in Action Potential Duration. A. Langendorff-perfusion set-up for MAP measurements. B. Representative MAP tracings of individual LQT2 (red) and wildtype (LMC, black) rabbits. C. Bull's eye plots of average AP D 75 (ms). Circles represent (from outside to inside) LV base, mid, and apex. APDs are color-coded as indicated. Septal APDs were not assessed. D. APD-Dispersion (SD of APD) in base, mid, apex, and global LV at 2, 3, and 4 Hz. E. ECG parameters in LQT2, LMC, and E4031-treated (1µg/kg/min) rabbits. n=10 LQT2, n=9 LMC, and n=10 E4031 rabbits; * p<0.05, * * p<0.01, * * * p<0.001 vs. LMC

  Figure B.1: Modeling the anatomy of patient Case #1 (from left to right): extracting geometry from MR datasets, generating volumetric meshes.

  Figure B.2: Estimation of the motion of Volunteer #1 (from left to right): end-diastolic frame, end-systolic frame.

  Fig.B.3 for a LBBB case (depolarization wave starting on the right ventricle) after registration and interpolation of the Ensite non contact mapping data available only on the left ventricle.

Figure B. 3 :

 3 Figure B.3: Personalized electrophysiology simulation obtained after registration and interpolation of the EnSite Depolarization Time isochrones available on the Left Ventricle. (For interpretation of the color map in this figure, the reader is referred to the web version of this article.)

  Figure B.5: Schematic representation of the Unscented Transform algorithm.

  Figure B.6: (Top) Bull's eye plots of the AHA 17 segment division (View from the apex).

  Figure B.7: Division of the myocardium into 6 regions for the left ventricle and 1 for the right ventricle, combining the standard AHA zones. (View from apex)

Figure B. 8 :

 8 Figure B.8: Comparison of the personalized motion at end-systole on short axis view fromthe base (left) and long axis (right) plane. Grey surfaces represent the ground truth, green contours the initial simulation, pink the estimated contours using regional volumes and blue the estimated contours using surface points. Since no clear difference can be made from the pink and blue contours, we conclude that both techniques give similar results and manage to match the ground truth. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

  Figure B.10: Comparison of the errors using the two types of observations, for various contractilities.

Figure B. 17 :

 17 Figure B.17: Case #1: Evolution of the estimated relative contractilities along the personalization steps. We note a mean relative contractility of 0.55. Contractility of zone #17 might be overestimated since it corresponds to the apex which is constrained by our boundary conditions.

  Fig.B.20.

Figure B. 19 :

 19 Figure B.19: Estimated global parameters for the 3 pathological cases compared to the healthy controls (box plots). σ 0 is the global contractility, µ the friction parameter, K the Bulk modulus related to the stiffness of the passive myocardium and R p the peripheral resistance in the aorta.

Figure B. 20 :

 20 Figure B.20: Registered outward flow curves for the pathological cases compared to the mean healthy curve, showing that Case 3 only has a smaller inflow than outflow.

Fig

  Fig. B.22 shows a rheological representation of this model.

Figure B. 22 :

 22 Figure B.22:Full electromechanical model. W e is the strain energy of the extracellular matrix considered here as an isotropic material, associated with a dissipative term η. u is a control variable which is driven by changes in transmembrane potential. It controls the contraction stress τ c . µ deals with the friction in the sarcomere while E s is a linear spring to enforce elasticity of the Z-discs (titin).

Figure C. 1 :

 1 Figure C.1: Examples and mutual positioning of the short and the main left ventricular long axis cardiac MR views used in this paper. See also C.8 which illustrates inter-subject variability of these views.

  Figure C.3: Classification forest pipeline for view recognition from image miniatures [Margeta 2014a]. Discriminative pixels from image miniatures are chosen from a random pool as features for a classification forest. The training dataset is augmented to improve robustness to the differences in the acquisitions without the need for extra data and labels.

Figure C. 4 :

 4 Figure C.4: Averages of the DETERMINE dataset for each cardiac view after cropping the central square. A reasonable alignment of cardiac structures can be observed in the dataset. The cardiac cavities and main thoracic structures can be seen.

  Figure C.5: Image miniatures by min-and max-pooling local intensities [Margeta 2014a].

Figure C. 6 :

 6 Figure C.6: Our CardioViewNet is based on CaffeNet network structure and is adapted for cardiac view recognition. We initialized the network weights from a pretrained CaffeNet. We then replaced the last 1000-dimensional logistic regression layer (previous fc8) with a 5-dimensional for 5 cardiac views. Then we finetuned the network with our data. We also extracted features from the last 4096-dimensional fully connected layer fc7 (in dark gray) from both CaffeNet and our finetuned Car-dioViewNet and used them with a linear SVM and a decision forest classifier. We achieve the best performance with our finetuned network directly. ReLU -rectified linear unit, max -max pooling, LRN -local response normalisation, conv-convolutional layer, fc -fully connected layer.

  Figure C.8: Typical examples of the training images from the DETERMINE dataset for different views. Note the acquisition and patient differences. In addition, the short axis slices cover the heart from the apex to the base with quite different appearances.

Figure C. 9 :

 9 Figure C.9: Sum of the confusion matrices over the 10 folds of our crossvalidation on the DETERMINE dataset for the best model classifying images using DICOM normal information and the best image-based predictor (using our finetuned neural network).

Figure

  Figure C.10: Examples for some of the least confident correct (the normal predictions are usually very peaky) predictions using CardioViewNet. Predicted and true (in parentheses) labels shown under the images. Below them are view-specific prediction confidences.

Figure C. 11 :

 11 Figure C.11: Example misclassifications using CardioViewNet. Predicted and true label (in parentheses) are indicated under the images and below them are viewspecific prediction confidences. The failures typically happen with less typical acquisitions and truly ambiguous images. The misclassified 2CH image indeed looks like a 3CH image except that the left atrium is missing and the ventricle is acquired at a non standard oblique angle.Similarly for the 4CH view, the right atrium is missing and one can already see parts of the outflow tract branching out of the left ventricle typical for a 3CH view. The 3CH view is captured with a detailed region of interest not very common in the dataset. Extra data augmentation could probably help to fix this case. Finally, the LVOT views are severely under-represented and can be confused with basal short axis views. Note that for all of these cases, the correct prediction is in the top 2 most likely views.

  To achieve great things, two things are needed: a plan and not quite enough time. LEONARD BERNSTEIN Clinical Context . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Manuscript Organization . . . . . . . . . . . . . . . . . . . . . 3 1.3 Acknowledgment . . . . . . . . . . . . . . . . . . . . . . . . . . 4

	i
	Chapter 1
	Introduction
	Contents
	1.1

Radiofrequency Ablation Planning for Cardiac Arrhythmia Treatment using Biophysical Modeling and Machine Learning Approaches

Modified from http://tapastic.com/episode/142659

Table 2 .

 2 1: Strategies for substrate-guided VT ablation.

	Ablation Lines
	Linear ablation from		Ablation Design
	Pacemapping	Isthmus Identification	Manoeuvers for
	Scar border zone	Ablation Target	Substrate
	Non-inducibility	Endpoing	Procedure
	+ + ++	of RF	Extension

Table 3 .

 3 1: Baseline patient characteristics and MRI findings.

	Patient	1	2	3	4	5	6	7
	Condition	ICM	ICM	ICM	ICM	ICM	NICM	NICM
	Gender	M	M	M	M	F	M	M
	Age	73	69	64	60	65	75	81
	Co-morbidities	HTN	None	HTN, AF	DM	None	AF	PAF
	B-blocker	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	ACE-I/ARB	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	Statins	Yes	Yes	Yes	Yes	Yes	No	Yes
	LVEF, %	27	35	35	25	31	36	17
	LVEDV, ml	199	292	245	304	185	196	285
	LV mass, g	129	172	147	182	118	89	206
	Scar on MRI	Yes	Yes	Yes	Yes	Yes	Yes	Yes
	Scar core, g	8	22	20	12	16	7	12
	Gray zone, g	12	26	4	19	7	8	5

ACE-I indicates angiotensin converting enzyme inhibitor; AF, atrial fibrillation; ARB, angiotensin receptor blocker; MRI, magnetic resonance imaging; EDV, enddiastolic volume; EF, ejection fraction; ICM, ischaemic cardiomyopathy; LV, left ventricle/ventricular; NICM, non-ischaemic dilated cardiomyopathy, PAF, paroxysmal atrial fibrillation.

Table 4 .

 4 1: Patient Electrophysiological Mapping Dataset Characteristics.

		# Points	# Points	# Points Map Source # Healthy # LAVA
		Original Post-Filtering On Scar		On Scar	On Scar
	P1	175	100	91	endocardium	54	37
	P2	700	522	83	endocardium	50	33
	P3	322	293	124	epicardium	113	11

Table 4 .

 4 2: Area under the curves for Precision-Recall and ROC.

		PR AUC
		Patient 1	Patient 2	Patient 3
	Full Features	0.79	0.88	0.75
	50% MRF	0.72	0.86	0.72
	25% MRF	0.73	0.91	0.81
	10 MRF	0.69	0.80	0.66
	Intensity Features	0.78	0.75	0.31
		ROC AUC
		Patient 1	Patient 2	Patient 3
	Full Features	0.80	0.89	0.92
	50% MRF	0.76	0.86	0.91
	25% MRF	0.73	0.91	0.91
	10 MRF	0.73	0.81	0.88
	Intensity Features	0.73	0.78	0.63

Table 4 .

 4 3: Classification performance scores using sample uncertainty weights.

	Confidence-based Training. Area Under the PR Curves
		Patient 1	Patient 2	Patient 3
	No Confidence Weighting	0.79	0.88	0.75
	Projection Distance Weighting	0.85	0.94	0.96
	Temporal Displacement Weighting	0.86	0.95	0.95
	Combined Uncertainty Weighting	0.84	0.94	0.94
	Confidence-based Training. Area Under the ROC Curves
		Patient 1	Patient 2	Patient 3
	No Confidence Weighting	0.80	0.89	0.92
	Projection Distance Weighting	0.85	0.94	0.99
	Temporal Displacement Weighting	0.87	0.95	0.99
	Combined Uncertainty Weighting	0.87	0.94	0.99

Table 4 .

 4 4: Random forests construction parameters when performing nested cross validation with precision optimization and temporal displacement weighting.

		Random Forest Construction Parameters
		PR Score	Trees	No. Features	No. Inner Folds
	Patient 1	0.86	25	60	5
	Patient 2	0.95	25	60	5
	Patient 3	0.95	10	20	5
	Table 4.5: Confusion Matrix For nested cross-validation with precision optimization
	and temporal displacement weighting.
	Patient	TN (healthy)	TP (lava)	FN	FP
	Patient 1	43		27	10	11
	Patient 2	46		26	7	4
	Patient 3	112		5	6	1

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 Clinical Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2.1 Electrophysiological Data . . . . . . . . . . . . . . . . . . . . 68 5.2.2 Imaging Data . . . . . . . . . . . . . . . . . . . . . . . . . . . [Cabrera-Lozoya 2015a] Cabrera-Lozoya, R., Berte, B., Cochet, H., Haissaguerre, M., Jais, P., Ayache, N., Sermesant, M. "Image-based Simulation of LAVA Intracardiac Electrograms." to be submitted. Chapter 5. Image-based Simulation of LAVA Intracardiac Electrograms

	Chapter 5
	Image-based Simulation of LAVA
	Intracardiac Electrograms
	Contents
	5.1

69 5.3 Cardiac Electrophysiology (EP) Simulation . . . . . . . . . . 69 5.3.1 Pseudo-personalisation: Zone-specific EP Properties . . . . . 70 5.3.2 Sinus Rhythm Computation from Surface ECGs . . . . . . . 71 5.4 Intracardiac EGM Simulation . . . . . . . . . . . . . . . . . . 72 5.4.1 The Dipole Approach . . . . . . . . . . . . . . . . . . . . . . 73 5.4.2 Unipolar and Bipolar EGM Computations . . . . . . . . . . . 75 5.5 Signal Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.6 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 77 5.6.1 CARTO Signal Characterization . . . . . . . . . . . . . . . . 77 5.6.2 Qualitative Assessment of Simulated EGMs . . . . . . . . . . 78 5.6.3 Simulated EGM Characterization . . . . . . . . . . . . . . . . 82 5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

Based on:

Table 5

 5

	.1: Simulation Parameter Values
	Parameter	Healthy Tissue	Grey Zone	Scar Tissue
	τ in [ms]	0.3	0.42	10×10 3
	τ out [ms]	6	8.4	10×10 3
	τ close [ms]	150	187.5	187.5
	τ open [ms]	120	120	120
	u gate	0.13	0.13	0.13
	Action Potential Peak Amplitude	1.0	0.7	0.7
	Conductivity	4.0	0.4	0.4

Table 5

 5

	.2: Structural Mesh Information and Electrode Locations per Patient
		P1	P2	P3	P4	P5
	# of Electrode Locations	1472	4804	1956	856	2800
	# of Tetrahedra in Mesh	112,308	96,260	150,480	96,288	29,586
	# of Nodes	23,355	20,306	31,100	21,231	6,742

Table 5

 5

		.3: KS test Results on Simulated Signals
	Feature	Patient 1 Patient 2 Patient 3 Patient 4 Patient 5
	Signal Range	0.009	<0.0001	0.02	0.002	0.01
	# Inflection Points	<0.0001	0.0003	<0.0001	<0.0001	<0.0001
	Signal Energy	0.04	-	<0.0001	-	<0.0001
	Dominant Frequency	0.01	<0.0001	0.005	-	-
	Mean Slope	<0.0001	0.002	<0.0001	0.0003	0.0005
	Fractionation Index	<0.0001	<0.0001	-	-	-
	Minimum Signal Value	-	0.0002	0.02	-	0.004
	Maximum Signal Value	0.004	<0.0001	0.0005	<0.0001	0.009

Figure 5.11: KS test Results on Simulated Signals.

Table 5 .

 5 4: EGM Characteristics of Clinical CARTO Signals

	Feature	Patient	non-LAVA	LAVA	p Value
			(µ ± σ)	(µ ± σ)	(CARTO)
	Signal Range	P1	5.31 ± 15.18	2.36 ± 1.73	<0.0001
		P2	1.40 ± 5.79	0.73 ± 0.42	<0.0001
		P3	3.84 ± 2.99	0.75 ± 0.42	<0.0001
		P4	2.09 ± 2.26	0.63 ± 0.54	<0.0001
		P5	3.79 ± 15.02	2.57 ± 0.48	0.001
	# Inflection Points	P1	66.97 ± 73.24	139.43 ± 104.456	<0.0001
		P2	315.47 ± 281.07	1085.6 ± 140.59	<0.0001
		P3	73.02 ± 102.99	246.02 ± 134.83	<0.0001
		P4	29.66 ± 90.22	258.22 ± 191.92	<0.0001
		P5	31.9 ± 43.97	51.74 ± 66.44	<0.0001
	Signal Energy	P1	(6.12 ± 4.39) ×10 -4 (6.68 ± 3.56) ×10 -4	0.02
		P2	(2.38 ± 5.14) ×10 -5 (8.90 ± 5.22) ×10 -4	<0.0001
		P3	(5.04 ± 3.59) ×10 -4 (6.67 ± 3.06) ×10 -4	<0.0001
		P4	(8.53 ± 4.48) ×10 -4 (6.65 ± 3.95) ×10 -4	0.002
		P5	(9.02 ± 2.50) ×10 -4 (9.27 ± 1.29) ×10 -4	<0.0001
	Dominant Frequency P1	17.09 ± 3.59	16.77 ± 3.61	-
		P2	20.46 ± 3.55	16.48 ± 3.65	<0.0001
		P3	16.73 ± 3.51	13.98 ± 4.12	<0.0001
		P4	17.04 ± 3.71	15.58 ± 4.6	-
		P5	17.65 ± 2.98	20.13 ± 2.87	<0.0001
	Mean Slope	P1	(4.19 ± 1.99) ×10 -3 (6.53 ± 2.88) ×10 -3	<0.0001
		P2	(1.21 ± 1.27) ×10 -2 (4.96 ± 0.74) ×10 -2	<0.0001
		P3	(3.92 ± 2.82) ×10 -3 (7.32 ± 2.44) ×10 -3	<0.0001
		P4	(4.80 ± 2.70) ×10 -3 (8.73 ± 3.80) ×10 -3	<0.0001
		P5	(3.44 ± 0.99) ×10 -3 (4.12 ± 1.75) ×10 -3	<0.0001
	Fractionation Index	P1	10.28 ± 6.42	16.21 ± 9.28	<0.0001
		P2	20.54 ± 14.18	24.07 ± 4.44	<0.0001
		P3	7.93 ± 5.39	16.46 ± 10.84	<0.0001
		P4	10.27 ± 6.76	13.84 ± 6.84	0.009
		P5	6.54 ± 2.53	8.54 ± 6.58	0.0003
	Minimum Value	P1	-2.82 ± 8.09	-1.06 ± 0.68	<0.0001
		P2	-0.68 ± 2.89	-0.38 ± 0.27	<0.0001
		P3	-1.65 ± 1.47	-0.36 ± 0.21	<0.0001
		P4	-1.08 ± 1.30	-0.34 ± 0.33	<0.0001
		P5	-1.67 ± 7.52	-1.07 ± 0.21	<0.0001
	Maximum Value	P1	2.49 ± 7.18	1.30 ± 1.18	0.01
		P2	0.72 ± 2.90	0.35 ± 0.25	<0.0001
		P3	2.18 ± 1.87	0.38 ± 0.26	<0.0001
		P4	1.01 ± 1.12	0.29 ± 0.24	<0.0001
		P5	2.11 ± 7.49	1.50 ± 0.28	0.007

Table 6 .

 6 1 and Figure6.7 summarize the statistics on the results of the different classification schemes across the five patients. Forest construction parameters are detailed in Table6.2.

Table 6 .

 6 For illustrative purposes, misclassified regions in Patient 1 after the fused-learning schemes were further analyzed. Because most of them correspond to false negatives areas (true LAVA which were classified as non-LAVA), we first assess the LAVAness of the CARTO clinical signals in these regions which were categorized as LAVA by an experienced electrophysiologist. The characteristics these electrograms were Biophysical Modeling Figure 6.7: Classification Statistics Summary On Five Patients.

		1: Classification Statistics On Five Patients
			Accuracy Sensitivity Specificity Precision
	Image-based	P1	94.0	73.3	98.5	91.6
		P2	92.6	33.3	100	100
		P3	91.8	39.0	99.0	90.0
		P4	93.7	11.0	100	100
		P5	92.5	10.0	100	100
		Mean	92.9	33.2	99.4	96.4
	Simulated EGM-based	P1	89.2	73.3	92.7	68.7
		P2	95.4	66.6	98.0	88.0
		P3	92.6	78.0	98.0	90.0
		P4	99.2	89.0	100	100
		P5	95.7	60.0	100	100
		Mean	94.4	73.2	97.8	89.4
	Maximum Confidence	P1	91.6	66.7	97.1	83.3
	Fusion	P2	92.6	33.0	100	100
		P3	96.7	78.0	99.0	94.0
		P4	97.6	67.0	100	100
		P5	96.8	70.0	100	100
		Mean	95.0	63.0	99.2	95.4
	Image-based,	P1	95.2	73.3	100	100
	Simulation-Enhanced	P2	92.6	33.3	100	100
		P3	96.2	74.0	99.0	94.0
		P4	97.6	67.0	100	100
		P5	93.6	40.0	100	100
		Mean	95.0	57.4	99.8	98.8
	Feature Set Fusion	P1	95.2	80.0	98.5	92.3
		P2	96.3	66.6	100	100
		P3	96.2	87.0	98.0	83.0
		P4	99.2	89.0	100	100
		P5	98.9	90.0	100	100
		Mean	97.2	82.4	99.2	95.0
	6.5.7 Ground Truth Re-evaluation

Table B . 1 :

 B1 Parameter values.

	Notation	Parameter Name	Value
	σ 0 (MPa)	Max Contraction	7
	k 0 (MPa)	Max Stiffness	5
	k AT P (s -1 )	Contraction Rate	15
	k RS (s -1 )	Relaxation Rate	35
	E s (MPa)	Linear Modulus	10
	α	Cross-bridges Unfasten Rate	0.5
	µ (MPa.s)	Viscosity	1.3
	c 1 (kPa)	Mooney Rivlin Modulus	10
	c 2 (kPa)	Mooney Rivlin Modulus	10
	K (MPa)	Bulk Modulus	10
	τ (s)	Wind. Charact. Time	0.8
	R p (MPa.m -3 .s) Wind. Periph. Resistance	100
	Z c (MPa.m -3 .s) Wind. Charact. Resistance	1
	L (kPa.s 2 .m -3 )	Wind. Total Art. Inertance

  . (Cited on page 111.) [Cherry 2004] Elizabeth M Cherry and Flavio H Fenton. Suppression of alternans and conduction blocks despite steep APD restitution: electrotonic, memory, and conduction velocity restitution effects. American Journal of Physiology-Heart and Circulatory Physiology, vol. 286, no. 6, pages H2332-H2341, 2004. (Cited on page 28.) [Chinchapatnam 2008] Phani Chinchapatnam, Kawal S Rhode, Matthew Ginks, C Aldo Rinaldi, Pier Lambiase, Reza Razavi, Simon Arridge and Maxime Sermesant. Model-based imaging of cardiac apparent conductivity and local conduction velocity for diagnosis and planning of therapy. Medical Imaging, IEEE Transactions on, vol. 27, no. 11, pages 1631-1642, 2008. (Cited on page 33.) [Chung 2008] E.S. Chung, A.R. Leon, L. Tavazzi, J.P. Sun, P. Nihoyannopoulos, J. Merlino, W.T. Abraham, S. Ghio, C. Leclercq and J.J. Bax. Results of the Predictors of Response to CRT (PROSPECT) trial. Circulation, vol. 117, no. 20, pages 2608-2616, 2008. (Cited on page 120.) [Ciaccio 2007] Edward J Ciaccio, Hiroshi Ashikaga, Riyaz A Kaba, Daniel Cervantes, Bruce Hopenfeld, Andrew L Wit, Nicholas S Peters, Elliot R McVeigh, Hasan Garan and James Coromilas. Model of reentrant ventricular tachycardia based on infarct border zone geometry predicts reentrant circuit features as determined by activation mapping. Heart Rhythm, vol. 4, no. 8, pages 1034-1045, 2007. (Cited on page 46.) [Ciaccio 2014] E.J. Ciaccio, H. Ashikaga, J. Coromilas, B. Hopenfeld, D.O. Cervantes, A.L. Wit, N.S. Peters, E.R. McVeigh and H. Garan. Model of bipolar electrogram fractionation and conduction block associated with activation wavefront direction at infarct border zone lateral isthmus boundaries. Circulation: Arrhythmia and Electrophysiology, 2014. (Cited on page 68.) [Ciresan 2012] D. Ciresan, U. Meier and J. Schmidhuber. Multi-column deep neural networks for image classification. In Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on, pages 3642-3649, June 2012. (Cited on page 148.) [Clancy 1999] Colleen E Clancy and Yoram Rudy. Linking a genetic defect to its cellular phenotype in a cardiac arrhythmia. Nature, vol. 400, no. 6744, pages 566-569, 1999. (Cited on page 28.) [Daelemans 2003] Walter Daelemans, Véronique Hoste, Fien De Meulder and Bart Naudts. Combined optimization of feature selection and algorithm parameters in machine learning of language. In Proceedings of the 14th European Conference on Machine Learning, pages 84-95. Springer, 2003. (Cited on page 21.) [D'Alonzo 1999] Albert J D'Alonzo, Jia L Zhu and Raymond B Darbenzio. Effects of class III antiarrhythmic agents in an in vitro rabbit model of spontaneous torsades de pointe. European Journal of Pharmacology, vol. 369, no. 1, pages 57-64, 1999. (Cited on page 107.) [De Bakker 1988] JM De Bakker, FJ Van Capelle, MJ Janse, AA Wilde, R Coronel, AE Becker, KP Dingemans, NM Van Hemel and RN Hauer. Reentry as a cause of ventricular tachycardia in patients with chronic ischemic heart disease: electrophysiologic and anatomic correlation. Circulation, vol. 77, no. 3, pages 589-606, 1988. (Cited on page 10.) [de Bakker 2010] J.M.T. de Bakker and F.H.M. Wittkampf. The pathophysiologic basis of fractionated and complex electrograms and the impact of recording techniques on their detection and interpretation. Circulation: Arrhythmia and Electrophysiology, 2010. (Cited on pages 68 and 72.) [De Chillou 2014] Christian De Chillou, Laurent Groben, Isabelle Magnin-Poull, Marius Andronache, Mohamed Magdi Abbas, Ning Zhang, Ahmed Abdelaal, Sonia Ammar, Jean-Marc Sellal and Jérôme Schwartz. Localizing the critical isthmus of postinfarct ventricular tachycardia: the value of pace-mapping during sinus rhythm. Heart Rhythm, vol. 11, no. 2, pages 175-181, 2014. (Cited on page 17.) [De Craene 2012a] Mathieu De Craene, Gemma Piella, Oscar Camara, Nicolas Duchateau, Etelvino Silva, Adelina Doltra, Jan D'hooge, Josep Brugada, Marta Sitges and Alejandro F Frangi. Temporal diffeomorphic free-form deformation: Application to motion and strain estimation from 3D echocardiography. Medical Image Analysis, vol. 16, no. 2, pages 427-450, Feb 2012. (Cited on pages 120 and 123.) [De Craene 2012b] Mathieu De Craene, Catalina Tobon-Gomez, Constantine Butakoff, Nicolas Duchateau, Gemma Piella, Kawal S Rhode and Alejandro F Frangi. Temporal Diffeomorphic Free Form Deformation (TDFFD) Applied to Motion and Deformation Quantification of Tagged MRI Sequences. In Statistical Atlases and Computational Modeling of the Heart (STACOM), volume 7085 of Lecture Notes in Computer Science, pages 68-77, 2012. (Cited on page 123.) [Decker 2010] Keith F Decker and Yoram Rudy. Ionic mechanisms of electrophysiological heterogeneity and conduction block in the infarct border zone. Ameri-

Verdandi is an opensource data assimilation library available at http://verdandi.gforge.inria.fr/. [Chapelle

2012a]

http://tetgen.berlios.de

http://www.hpfem.jku.at/netgen

SOFA is an open source medical simulation software (available at www.sofa-framework.org)

Acknowledgments Acknowledgments

We used data and infrastructure made available through the Cardiac Atlas Project (www.cardiacatlas.org - [Fonseca 2011]). See [Kadish 2009], [Tobon-Gomez 2013] for more details on the datasets.

DETERMINE was supported by St Jude Medical, Inc; and the National Heart, Lung and Blood Institute (R01HL91069). A list of participating DETERMINE investigators can be found at http://www.clinicaltrials.gov. This work uses scikit-learn toolkit ([Pedregosa 2011]) for decision forests and Caffe deep learning framework ( [Jia 2014]) for training of the convolutional neural network and the pretrained model (CaffeNet). This model was trained on a subset [Russakovsky 2014] of the ImageNet [Deng 2009] dataset.

funded by the European Research Council through the ERC Advanced Grant MedYMA 2011-291080 (on Biophysical Modeling and Analysis of Dynamic Medical Images).

Funding

This work was supported by Microsoft Research through its PhD Scholarship Programme and ERC Advanced Grant MedYMA 2011-291080. The research leading to these results has received funding from the European Union's Seventh Framework Programme for research, technological development and demonstration under grant agreement no. 611823 (VP2HF).

Based on a collaboration with researchers in the Department of Cardiology and Angiology I at the Heart Center Freiburg University, where I was responsible for the design and implementation of the machine learning approach. This work was accepted in the Hearth Rhythm Journal.

Odening, Katja E., et al. "Spatial correlation of action potential duration and diastolic dysfunction in transgenic and drug-induced LQT2 rabbits." Heart Rhythm 10.10 (2013): 1533-1541. Transgenic and Drug-induced LQT2 Rabbits none of the clinical manifestations of LQTS has been directly, patho-physiologically linked to the presence of severe mechanical dysfunction. Importantly, however, we could demonstrate that mechanical cardiac function was sufficiently altered in LQT2 to allow a differentiation between LQT2 and normal hearts solely based on phase contrast MRI. These data indicate that LQTS should probably not only be considered an electrical but rather an electromechanical disorder.

A.5 Conclusions

Prolongation and increased dispersion of cardiac repolarization led to globally and regionally impaired diastolic and systolic mechanical function and increased mechanical dispersion in transgenic and drug-induced LQT2 rabbits. Moreover, in transgenic LQT2 rabbits, regional differences in APD correlated with the extent of reduction of regional diastolic peak velocities and regional time-to-diastolic peak velocitiesas marker of impaired diastolic function, indicating that LQTS is not purely an electrical but rather an electromechanical disorder. Finally, due to this association of electrical and mechanical functions, it was possible to (1) differentiate between LQT2 and normal LMC rabbits' hearts, and to (2) discern differences in the extent of APD prolongation and APD dispersion in LQT2 rabbits solely based on MRI data, indicating potential new approaches for risk stratification and treatment decisions of LQTS patients.

A.6 Supplementary data

Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.hrthm.2013.07.038.

Appendix A. Spatial Correlation of APD and Diastolic Dysfunction in

Transgenic and Drug-induced LQT2 Rabbits

Appendix B

Personalization of a Cardiac Electromechanical Model using Reduced Order UKF from Regional Volumes

Based on a collaboration with another PhD student in the Asclepios Research Team, where I performed the electrophysiology personalization. This work was accepted in the Medical Image Analysis Journal.

Marchesseau, Stéphanie, et al. "Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes." Medical image analysis 17.7 (2013): 816-829. using Reduced Order UKF from Regional Volumes

B.1 Abstract

Patient-specific cardiac modelling can help in understanding pathophysiology and therapy planning. However it requires to combine functional and anatomical data in order to build accurate models and to personalize the model geometry, kinematics, electrophysiology and mechanics. Personalizing the electromechanical coupling from medical images is a challenging task. We use the Bestel-Clément-Sorine (BCS) electromechanical model of the heart, which provides reasonable accuracy with a reasonable number of parameters (14 for each ventricle) compared to the available clinical data at the organ level. We propose a personalization strategy from cine MRI data in two steps. We first estimate global parameters with an automatic calibration algorithm based on the Unscented Transform which allows to initialize the parameters while matching the volume and pressure curves. In a second step we locally personalize the contractilities of all AHA (American Heart Association) zones of the left ventricle using the Reduced Order Unscented Kalman Filtering on Regional Volumes. This personalization strategy was validated synthetically and tested successfully on eight healthy and three pathological cases.

B.2 Introduction

Understanding pathologies, planning a therapy, training for a surgery or selecting suitable patients for a treatment are challenges that cardiologists face everyday. For instance, Cardiac Resynchronization Therapy (CRT), which releases heart failure with a pacemaker, is still difficult to predict in spite of many advances [Chung 2008, Cazeau 2001], and fails in 30% of the treated cases. To ease their tasks, a significant amount of image data is now available. One goal of computational cardiology is to use cardiac imaging data to build patient-specific electromechanical models of the heart that can then be used to better analyze and simulate cardiac function [Smith 2011[START_REF] Sermesant | [END_REF].

In order to reach such a goal, several components are mandatory. Firstly, an accurate anatomical model has to be extracted from the image data to represent the heart geometry. Recently, various approaches have been proposed to automatically segment cardiac images from different modalities, such as echocardiography [Butakoff 2011[START_REF] Zhu | [END_REF], computed tomography (CT) [Ecabert 2011, Zheng 2008] and magnetic resonance imaging (MRI) [START_REF] Petitjean | [END_REF], Schaerer 2010]. In particular, MRI has the advantage of providing a good image quality noninvasively and without radiation. In this paper, we extract personalized geometrical models from MRI by using a 3D Active Shape Model approach by [Tobon-Gomez 2012]. The estimated motion of the myocardium is also needed to build patient-specific simulations. This motion can be estimated from sequences of MRI data using various algorithms [McLeod 2012[START_REF] Mansi | [END_REF], De Craene 2012a]. Then, an electrophysiological model has to be considered in order to simulate the electrical wave propagation within the myocardium which leads afterwards to the contraction of using Reduced Order UKF from Regional Volumes several passive and active terms, as well as some boundary constraints that represent the fluid flow that circulates inside the ventricles during the four cardiac phases [Marchesseau 2012a]. This model depends on 14 global parameters that need to be estimated. To this end, we propose to first assess global values for the parameters using the Unscented Transform algorithm [Julier 1997] from global observed quantities (endocardial volume and pressure) as described in [Marchesseau 2012b]; and then to use the ROUKF (Reduced Order Unscented Kalman Filter) personalization algorithm [Moireau 2011] implemented in Verdandi [Chapelle 2012a] from local measurements to assess more precisely the contractility per region (see Fig. C.3 for a representation of this pipeline).

B.3.2.1 Calibration based on the Unscented Transform

The initialization of the parameters is a crucial prerequisite for the personalization algorithm to converge towards a relevant solution. In [Marchesseau 2012a], we showed that this can be performed in a fast and efficient way using the Unscented Transform algorithm to match indices on the volume and pressure curves (if available) leading therefore to the calibration of 4 to 7 parameters. This algorithm, represented in Fig. B.5, builds a covariance matrix Cov(θ, Z) between the relevant parameters which are spread around some initial parameter set θ 0 and the observations Z -in our case the minimum of the LV (Left Ventricle) volume and the minimum and the maximum of its derivative. The new parameters θ new are then found to minimize the difference between the mean simulated observations Z and the measured observations Z obs with

In one simulation, this algorithm enables to calibrate the main parameters chosen from a sensitivity analysis on the volume curve (and the pressure curve when using Reduced Order UKF from Regional Volumes results in cardiac personalization: [START_REF] Billet | Assimilation de données images pour la personnalisation d'un modèle électromécanique du coeur[END_REF], Delingette 2012[START_REF] Sainte-Marie | Modeling and estimation of the cardiac electromechanical activity[END_REF] that use the variational approach to estimate the contractility parameters or [Xi 2011, Chabiniok 2011] that preferred the sequential approach). We use in this paper the Unscented Kalman Filter in its reduced form, derived in [Moireau 2011] since we only estimate the parameters.

Unscented Kalman Filtering

The Unscented Kalman Filter can be seen as an iterative Unscented Transform described in B.3.2.1. It is based on three steps at each discrete time n, for the state vector X which consists in the state values (positions, velocities, ...) and the parameters θ.

• Sampling: Creates the sigma-points (minimal set of sample points around the mean) given an interpolation rule leading to I, using a state error covariance matrix P:

were X (i) is the ith sigma-points associated with vector X.

• Prediction: Estimates the next state as the mean of the states simulated with the sigma-points, and the next covariance as the covariance of the states simulated with the sigma-points:

where the mean E α and the covariance Cov α are calculated with the weights α i of the sigma-points, and noting X * the concatenation of the vectors X (i) .

• Correction: Updates the state vector and its covariance given the new sigmapoints and corresponding observations using the new observation operator H and the covariance of the noise W through the following calculations:

Reduced Order Unscented Kalman Filtering We aim, in this paper, at estimating only the p parameters of the model (in our case the 17 regional contractilities from the 17 AHA zones of the left ventricle). Therefore, we only generate 2p + 1 sigma-points and the covariance P is reduced to using Reduced Order UKF from Regional Volumes close to 1 for the scar zone means that the detection fails (which is the case for the zones with sizes 300 tetrahedra, 650 tetrahedra and 1300 tetrahedra). On the other hand, the zone of size 2000 tetrahedra has an estimated relative contractility much smaller than 1 which means that the detection succeeds. Therefore, the proposed personalization method based on regional volumes is able to detect zones with lower contractility as long as these zones are large enough (here 3% of the myocardium).

B.4.3 Real Cases Personalization

Personalizing real cases requires three steps:

• Selecting the region that has the smallest displacement (basal area or apex) in order to constraint this region as a boundary condition.

• Performing a calibration of 4 parameters (including contractility) for the left ventricle, based on the volume curve.

• Running the ROUKF personalization algorithm using regional volumes on the LV endocardium in order to estimate the 17 regional contractilities of the LV.

This pipeline was applied on 3 pathological cases and 8 healthy volunteers and enabled to draw preliminary conclusions on pathological cases versus healthy controls. Results of the full personalization strategy are detailed here for the pathological Case #1, then error measurements and estimated contractilities are given on all cases. The observation vector Z is of size the number of regional volumes and X contains all the state variables and the parameters to estimate (node positions, velocities, ...). The full pipeline takes about 2p + 2 = 36 times a simulation time on a single CPU (one simulation lasts about 10 minutes for 1 heart beat on a 60 000 tetrahedral mesh) where p is the number of estimated parameters (here 17 regional contractilities). However, full parallelization can be done leading to a total computation time of twice a simulation time (20 minutes). .17 gives the estimated relative contractilities along the personalization steps and shows that the estimation converges before the end of the personalization.

B.4.3.2 Error Measurements

In order to evaluate the accuracy of the personalized motion, we computed the errors previously defined on all the tested cases, before and after personalization. Fig. B.18 presents in (a) the resulting errors on the regional volumes εRegV ol which is what is actually minimized during the personalization, as well as its value at end-systole (b). Finally, as an indicator of the global motion on all the surfaces, the distance error on the contours εC is given in (c). using Reduced Order UKF from Regional Volumes

C.1 Abstract

In this paper we propose a convolutional neural network-based method to automatically retrieve missing or noisy cardiac acquisition plane information from magnetic resonance imaging (MRI) and predict the five most common cardiac views. We finetune a convolutional neural network (CNN) initially trained on a large natural image recognition dataset (Imagenet ILSVRC2012) and transfer the learnt feature representations to cardiac view recognition. We contrast this approach with a previously introduced method using classification forests and an augmented set of image miniatures, with prediction using off the shelf CNN features, and with CNNs learnt from scratch.

We validate this algorithm on two different cardiac studies with 200 patients and 15 healthy volunteers respectively. We show that there is value in finetuning a model trained for natural images to transfer it to medical images.

Our approach achieves an average F1 score of 97.66% and significantly improves the state of the art of image-based cardiac view recognition. This is an important building block to organise and filter large collections of cardiac data prior to further analysis. It allows us to merge studies from multiple centers, to perform smarter image filtering, to select the most appropriate image processing algorithm, and to enhance visualisation of cardiac datasets in content based image retrieval. Acquisition Plane Recognition Figure C.2: DICOM plane normals for different cardiac views [Margeta 2014a]. In our dataset distinct clusters can be observed (best viewed in color). Nevertheless the separation might not be the case for a more diverse collection of images. Moreover as we cannot always rely on the presence of this tag, so an image-content based recognizer is necessary.

the second group learn to predict the target labels from the image content (pixel intensities) directly without the need for DICOM tags. This group includes the random forest classifier with image miniatures (C.3.2) and all convolutional neural network approaches (C.3.3). Apart from the view label linked to each image in the training set, no other information is necessary to train the classification models. To increase the number of the training samples (for image content based algorithms) we augment the dataset with small label preserving transformations such as image translations, rotations, and scale changes. See C.3.2.2 for more details.

C.3.1 Using DICOM Meta-information: Plane normal + Forest

Both [Zhou 2012] and [Margeta 2014a] showed that where the DICOM orientation (0020,0037) tag is present we can use it to predict the cardiac image acquisition plane (See Fig. C.2). This tag is not defined as a cardiac view but as two 3dimensional vectors defining orientation of the imaging plane with respect to the MR scanner coordinate frame. It is straightforward to compute the 3-dimensional normal vectors of this plane as a cross-product of these two vectors specified in the tag. We then feed these three-dimensional feature vectors into any classifier, in our case a classification forest [START_REF] Breiman | Random forests-random features[END_REF]]. This method is shown in the results section as Plane normal + forest.

This method uses feature vectors computed from the DICOM orientation tag and cannot be used in the absence of this tag. This happens for example after Acquisition Plane Recognition output of a neuron to be 0.5. The finetuning is quite efficient and takes approximately 4 hours on a single NVIDIA Tesla M2050 GPU for 8000 iterations. Results of this method are presented as CardioViewNet. We also show results for prediction of an SVM classifier using fc7 features extracted from the finetuned network as CardioViewNet fc7 + SVM. In other words, we replace the 1000-class multinomial regression layer of the finetuned network by a linear SVM classifier. Results using a classification forest instead are listed as CardioViewNet fc7 + forest. The possibility to replace the final classifier is important to quickly train new classifiers for additional views without extra finetuning. In addition to the training set augmentation, we perform oversampling at the test time i.e. average predictions of ten 227x227 image crops: the central crop and the four 227x227 corner aligned crops and their horizontally flipped versions (vertically flipped images are rare). We report these results as CardioViewNet + oversample. We will see that this can improve performance on an independent dataset.

C.3.3.3 Training the network from scratch

Good initialisation of the network is important and the CaffeNet trained on the ImageNet dataset helps us to get well started. The initial motivation behind the finetuning was that there were too few images in our dataset to train the whole network from scratch. While the number of images is certainly smaller than in the ImageNet dataset, our target domain is also much simpler. We are predicting only 5 classes whose appearance variability is lower than the one across the ImageNet classes (e.g. variability of felines in different poses and appearances when they are all labeled as a cat). To test whether there is any value in the finetuning instead of learning the network parameters from scratch we train from the ground up two networks. First, a simpler LeNet-5 network [START_REF] Lecun | [END_REF]] (shown as LeNet-5 from scratch) as defined by Caffe but with the last layer adapted to a 5 class Acquisition Plane Recognition DETERMINE KCL F1 score precision recall F1 score precision recall Acquisition Plane Recognition handcrafting the features or training with extra annotations. This also allowed us to gain performance over models learnt from scratch. However, even the performance of models learnt from scratch is very encouraging for further exploration.

Features extracted from our network should be useful as descriptors for new views and extend our method even to other pathology specific views (such as those used in congenital heart diseases) and acquisition sequences other than SSFP, and even to recognize the acquisition sequences themselves. The methods presented in this paper are important additions to the arsenal of tools for handling noisy metadata in our datasets and are already helping us to organize collections of cardiac images.

In the future this method will be used for semantic image retrieval and parsing of medical literature.