E Pernod

M Sermesant
email: maxime.sermesant@inria.fr

E Konukoglu

J Relan

H Delingette

N Ayache

A Multi-Front Eikonal Model of Cardiac Electrophysiology for Interactive Simulation of Radio-Frequency Ablation

Keywords: Medical interactive simulation, Modeling of the heart, Haptic device

Virtual reality based therapy simulation meets a growing interest from the medical community due to its potential impact for the training of medical residents and the planning of therapies. However, computer models of the human anatomy are often very computationally demanding, thus incompatible with the constraints of such interactive simulations.

In this paper, we propose a fast model of the cardiac electrophysiology based on an eikonal formulation implemented with an anisotropic fast-marching method. We demonstrate the use of this model in the context of a simulator of radio-frequency ablation of cardiac arrhythmia from patient-specific medical imaging data. Indeed, this therapy can be very effective for patients but still suffers from a rather low success rate. Being able to test different ablation strategies on a patient-specific model can have a great clinical impact.

In our setting, thanks to a haptic 3D user interface, the user can interactively measure the local extracellular potential, pace locally the myocardium or simulate the burning of cardiac tissue as done in radio-frequency ablation interventions.

Introduction

Cardiac arrhythmias are pathologies characterised by an abnormal heart beat or rhythm. The contraction and relaxation of the heart is indeed controlled by the depolarisation and repolarisation of the cardiac cells that propagates through the myocardium, from the sinus node, the natural pacemaker of the heart, to the ventricles. The abnormal propagation of the electric wave inside the heart leads to dysfunction of cardiac mechanics causing a loss of efficiency when pumping the blood to the general circulation.

Treatment of cardiac arrhythmias has considerably changed, for instance with radio-frequency (RF) ab-lation becoming widely available as an alternative to drug therapy. These interventions are carried out under x-ray guidance, with specialised catheters for making invasive recordings of the electrical activity in the heart, and even reconstruct the chamber geometry. RF ablation can be used to burn the pathological cardiac cells responsible for abnormal electrophysiological propagations. One source of difficulty for cardiologists when performing this intervention is to locate those cells based on the reading of multiple electrical signals acquired in the heart. In many cases, the therapist adopts a trial and error approach by continuously checking if the electric propagation is altered by the performed ablation. Furthermore, partly due to the challenge in performing the intervention, a large proportion (around 50 %) of patients being ablated have a recurrence of the pathology at least 6 months after the intervention.

Therefore, these procedures can be highly effective with minimal side effects but some groups of patients have unsatisfactory success rates and may entail long procedures with high x-ray radiation dose to both patient and staff. Moreover, serious side effects can arise if the lesions extend beyond the target area, for the ablations. There is a need for substantial innovation in order to reliably achieve successful results in an acceptable time, with lower radiation dose and reduced risk of accidental damage to adjacent structures.

The aim of this research work is to design models of the cardiac electrophysiology that are suitable for clinical use and to propose methods to combine these models with interventional data in order to better estimate the patient cardiac function and help in the training, planning and guidance of such procedures.

During the last decades, important progress has been made in modelling the cardiac electrophysiology, with the design of more and more complete ionic models [START_REF] Tusscher | Panfilov, A model of the human ventricular myocyte[END_REF]. However, the integration of such models at the scale of the organ often proves difficult, as this results in complex reaction-diffusion equations which require careful mathematical analysis and produce important computational load. See for instance [START_REF] Garny | Dimensionality in cardiac modelling[END_REF] for an analysis of this problem.

To tackle this, rule-based models were developed in many areas of computational biology allowing to study the spatial interaction of cells with discretised behaviours. This leads to a better understanding of interaction phenomena at a much larger scale. There have been many examples of cellular automata used to model excitable media. However, cellular automata do not derive from the discretisation of a continuous model, thus the impact of the topology and size of the grid on the resulting dynamics is difficult to control, for instance on the anisotropy [START_REF] Schnfisch | Anisotropy in cellular automata[END_REF]. Different extensions to the classical cellular automata have been proposed in the literature to better reproduce the behaviour of reaction-diffusion models (like curvature [START_REF] Gerhardt | A cellular automaton model of excitable media: II. curvature, dispersion, rotating waves and meandering waves[END_REF]) but these often remain ad hoc as this is not the discretisation of a continuous model.

In this article, we propose to integrate both the continuous dynamics, through an eikonal approximation of the reaction-diffusion systems, and the discrete states, through a cellular automaton, in order to achieve fast simulations of the propagation, while preserving most of the properties of the reaction-diffusion model. To achieve this, we use the Fast-Marching Method (FMM) with several improvements, in order to be able to take into account the anisotropy of the propagation and the existence of several fronts.

Virtual-reality based surgery simulation opens up new fields in medical training and planning of therapies. In this paper, we investigate the development of a simulator prototype for rehearsing RF ablation of the myocardium based on the proposed model.

Electrophysiology Models

Modelling the cell electrophysiology is an active research area since the seminal work of Hodgkin and Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. The precise modelling of the myocardium involves a cell membrane model embedded into a set of partial differential equations (PDE) modelling a continuum. We can divide these models into three categories, from the more complex to the simpler (numerically):

• Biophysical models: semi-linear evolution Partial Differential Equation (PDE) with ionic models (up to fifty equations for ions and channels [START_REF] Noble | A modification of the Hodgkin-Huxley equations applicable to purkinje fibre action and pace-maker potentials[END_REF][START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibers[END_REF][START_REF] Luo | A model of the ventricular cardiac action potential: depolarization, repolarization, and their interaction[END_REF][START_REF] Noble | Improved guinea-pig ventricular cell model incorporating a diadic space, I Kr and I Ks , and length and tension dependent processes[END_REF], 1])

• Phenomenological models: semi-linear evolution PDE with mathematical simplification of biophysical models (bidomain, monodomain [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF])

• Eikonal models: one static non-linear PDE for the depolarisation time derived from the previous models (Eikonal-Curvature [12], Eikonal-Diffusion [START_REF] Colli Franzone | Wavefront propagation in activation model of the anisotropic cardiac tissue: Asymptotic analysis and numerical simulations[END_REF])

PDE based electrophysiology models are more sophisticated and can account for complex physiological states such as ischemia or arrythmia. However, they often rely on a large set of parameters that are difficult to personalize and are computationally very demanding at the organ scale. Indeed the space scale of the electrical propagation fronts is in general much smaller than the size of the ventricles and fine mesh must be used.

On the other hand, eikonal models only simulate the time at which a front reaches a given point and therefore lead to much faster computations. Such models can be used to estimate patient-specific parameters [START_REF] Chinchapatnam | Model-based imaging of cardiac apparent conductivity and local conduction velocity for diagnosis and planning of therapy[END_REF], or explore a large set of simulations [START_REF] Camara | Personalization of fast conduction purkinje system in eikonal-based electrophysiological models with optical mapping data[END_REF].

Eikonal models are derived as asymptotic solutions of PDEs and therefore both approaches can be easily compared [START_REF] Jacquemet | An eikonal approach for the initiation of reentrant cardiac propagation in reaction-diffusion models[END_REF][START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF]. Differences can appear due to the boundary conditions: for PDEs, these are Neumann boundary conditions, which impose the no-flux condition at the boundary, whereas the eikonal equation has no such constraints. Another difference can arise from the fact that in PDEs the speed of the travelling front is time dependent whereas in eikonal form, only a constant speed asymptotic value is considered. In [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF], Konukoglu et al. have added correcting terms in the eikonal models in order to minimize those discrepancies.

One important limitation of eikonal models is that they only describe the time of depolarization without describing the state of repolarization and refractoriness of cardiac cells. Those phenomena are very important for the simulation of arrythmias. Jacquemet [START_REF] Jacquemet | An eikonal approach for the initiation of reentrant cardiac propagation in reaction-diffusion models[END_REF] proposed an extension of those models to consider periodic arrythmia without taking into account the dynamic aspect of restitution curves. In the remainder, we detail a multi-front eikonal model which captures the main characteristics of cardiac electrophysiology propagation (anisotropy, refractory period, repolarisation, restitution curve) at a low computational cost [START_REF] Sermesant | An anisotropic multifront fast marching method for real-time simu-lation of cardiac electrophysiology, in: Proceedings of Functional Imaging and Modeling of the Heart[END_REF].

The Anisotropic Fast Marching Method

The classical FMM [START_REF] Sethian | Level set methods and fast marching methods: Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF] can be used to very efficiently solve the standard isotropic eikonal equation:

F (X) ∇T (X) t D(X)∇T (X) = 1 (1)
where T (X) is the arrival time at point X and F (X) the local speed function. D(X) is a 3 × 3 symmetric tensor defining the anisotropy of the propagation at point X. In the classical isotropic case, D(X) is the identity matrix. Whereas in the anisotropic case, it captures the fact that the electrical propagation is roughly three times faster along the fibre orientation than in its transverse direction. Therefore, this tensor is computed as:

D(X) = 1 9 I + 8 9 f (X) ⊗ f (X) ( 2
)
where I is the identity matrix and f (X) is the unit fibre orientation vector at point X.

In this method, points are separated into three sets: U N KN OW N , T RIAL, and KN OW N , and the T RIAL set of points X is ordered along the increasing associated times T (X). N (X) is the neighbourhood of X. In the initialisation step, all points are labelled as U N KN OW N with a time of +∞, starting points are labelled as KN OW N with a time of 0, and their neighbours are iteratively updated. The FMM algorithm is described here briefly, see Algorithm 1.

Algorithm 1 Fast Marching Method

while TRIAL list is not empty do X ← argmin X∈T RIAL {T (X)} remove X from T RIAL and add X to KN OW N for all X i ∈ N (X) and

X i / ∈ KN OW N do compute T (X i ) = UPDATE(X i , X) if X i /
∈ T RIAL then add X i to T RIAL end for end while Taking into account the anisotropy of the propagation is important in many cases of excitable media. But introducing this phenomenon in the FMM is not trivial. Indeed, the regular FMM uses the collinearity between the time gradient and the characteristic direction to solve the equation very efficiently, but in the anisotropic case these directions do not match any more.

There have been many different algorithms proposed to solve such anisotropic equations using single-pass [START_REF] Sethian | Ordered upwind methods for static hamilton-jacobi equations: theory and algorithms[END_REF][START_REF] Kevorkian | Partial differential equations: Analytical solution techniques[END_REF] or iterative [START_REF] Kao | Fast sweeping methods for static hamilton-jacobi equations[END_REF][START_REF] Qian | A fast sweeping method for static convex hamilton-jacobi equations[END_REF] methods. We proposed a fast algorithm for solving anisotropic eikonal equation on general meshes without increasing the neighbourhood size around each vertex and following the characteristic direction similar to singlepass methods. The idea is to include a recursive correction scheme taking into account the fact that, due to anisotropy, the immediate neighbourhood used for computation may not always contain the characteristic direction at the time it is computed. To achieve this, an additional CHAN GED list is introduced which is also empty at the beginning and which is used during the recursive correction step. This efficient algorithm can cope with very important anisotropies and can be applied to more general forms of static convex Hamilton-Jacobi equations, and on Cartesian or unstructured grids [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF]. Pseudocode is provided in the following algorithms 2 and 4 (U P DAT E function in Annex).

The anisotropic FMM algorithm is not significantly longer to compute than the standard FMM algorithm, especially when the anisotropy is moderate. For instance with a factor of three between longitudinal and transverse speeds as in cardiac tissue, simulation on a 13 000 nodes mesh only takes few seconds. The results are presented on a cubic mesh, composed of 13 000 nodes and 75 000 tetrahedra. The initialisation is in the front top corner (see Figure 1).

A Multi-Front Fast Marching Approach

The FMM is a static method to solve for an evolving front, in the sense that the variable of the eikonal equation is time. But for our purpose, we can have several fronts evolving at the same time, because cardiac cells go back to a resting state, so they can be excited again before the first front disappeared. The key idea is to introduce a time-stepping scheme while using the FMM to compute the propagation during each time step.

In order to achieve this, we also have to introduce a refractory state, between the excited and the equilibrium state, otherwise, any node going back to the equilibrium would immediately be excited again by his neighbours still excited. This is in agreement with cardiac electrical action potential propagation, where cells have this refractory period [START_REF] Tusscher | Panfilov, A model of the human ventricular myocyte[END_REF]. The Multi-Front FMM is described in algorithm 3.

Algorithm 2 Anisotropic Fast Marching

while T RIAL or CHAN GED are not empty do if CHAN GED is not empty then

X ← argmin X∈CHAN GED {T (X)} remove X from CHAN GED else X ← argmin X∈T RIAL {T (X)} remove X from T RIAL and add X to KN OW N end if for all X i ∈ N (X) and X i ∈ KN OW N do compute T (X i ) = UPDATE(X i , X) if T (X i ) < T (X i ) then T (X i ) ← T (X i ) add X i to CHAN GED list end if end for for all Y i ∈ N (X) and Y i ∈ T RIAL ∪ U N KN OW N do compute T (Y i ) = UPDATE(Y i , X) if Y i ∈ T RIAL and T (Y i ) < T (Y i ) then T (Y i ) ← T (Y i ) else if Y i ∈ U N KN OW N then T (Y i ) ← T (Y i ) remove Y i from U N KN OW N and add Y i to T RIAL end if end for end while
This Multi-Front FMM is able to simulate an impulse train on a piece of excitable media (see Figure 2).

Cardiac Electrophysiology Simulation

The presented algorithm can be used to simulate transmembrane potential propagation in cardiac tissue. The idea is to compute the depolarisation time using the Multi-Front Anisotropic FMM. The local speed F corresponds to the local tissue conductivity. Then, when a given vertex is depolarised, we use the time-stepping to know when the vertex goes into repolarisation and then refractory period, and when it is excitable again. The time during which a cell is depolarised is called the Action Potential Duration (APD) and it corresponds to the period during which the cardiac fibres are actively mechanically stretched. Furthermore, we can compute a pseudo transmembrane potential, which is at 10 mV when the vertex is depolarised (red), and -90 mV when it is excitable (blue), and interpolated using a precomputed curve in between (see Figure 3). Compared to similar approaches in propagation simulation [START_REF] Reumann | Preventive ablation strategies in a biophysical model of atrial fibrillation based on realistic anatomical data[END_REF], it could be seen as a "continuous cellular automaton": the states are still discrete, but the propagation is continuous.

To better simulate the cardiac cells' behaviour, it Algorithm 3 Multi-Front Fast Marching Algorithm on a period T . IT : Integrated Time, ET : Elapsed Time, RT : Repolarisation Time, RP : Refractory Period IT = 0.0 while IT < T do ET = 0.0 while T RIAL not empty and ET < timeStep do X ← argmin X∈T RIAL {T (X)} remove X from T RIAL and add X to KN OW N for all X i ∈ N (X) and

X i ∈ U N KN OW N do compute T (X i ) = UPDATE(X i , X) if X i / ∈ T RIAL then add X i to T RIAL ET = T (X i ) -IT end for end while IT + = timeStep for all X ∈ KN OW N do if IT -T (X) > AP D(X) then remove X from KN OW N and add X to REF RACT ORY RT (X) ← T (X) + AP D(X) end if end for for all X ∈ REF RACT ORY do if IT -RT (X) > RP (X) then remove X from REF RACT ORY and add X to U N KN OW N end if end for end while
is important to introduce the restitution curve, which is the relationship between the pacing frequency and the APD [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF]. Several different curves can be found in the literature, depending mostly on the experimental protocol [START_REF] Nash | Whole heart action potential duration restitution properties in cardiac patients: a combined clinical and modelling study[END_REF]. This can be directly used within our model, as we compute the APD at cycle n + 1 as a function of the previous diastolic interval (DI) DI n (time during which a cell is repolarised):

AP D n+1 = T close * log 1 - (1 -h min ) * exp -DIn Topen h min (3)

Radio-Frequency Ablation Model

There are numerous publications on the detailed simulation of the radio-frequency ablation, as this techniques is widely used in different clinical applications (liver tumours, cardiac arrhythmias) [27, 28,

!"# $% #& $% !"# $'(% )" $%
)" $'(% !*+,-./0%

)/12-*.324%

)/5.% 29]. Such models aim at predicting the evolution of the tissue temperature as a function of the delivered energy taking into account various cooling effects. This is important in order to ensure that the heated tissue creating lesions do cover the targeted volume, in an optimal time.

Here, as a first approximation, we used a very simplified model for which the lesion size is user-defined and the conductivity F of tetrahedral elements considered as ablated is set to zero.

The parameter values of the electrophysiology model as function of tissue type are summarized in table 1. Restitution curves for healthy and peri-infarct cardiac tissues are shown in Figure 4.

Implementation

Visualisation and Computation in SOFA

The proposed work has been implemented in an open source framework, SOFA [START_REF] Allard | Sofa an open source framework for medical simulation[END_REF] suitable for realtime modelling of deformable structures and mechanical behaviours, particularly for medical simulation. SOFA is an open-source library mostly written in C++ and using XML script language. Its architecture has been designed to be very modular with a lean core containing the overall architecture and prototype functions and on the other side, a large set of modules derived from the core. To link the different components (C++ modules of the framework) together, scene graphs (simulation tree in XML format) are used to describe a simulation. In this work, several modules related to haptic rendering and collision detection have been used to provide a satisfactory user interface. Furthermore, we developed specific data containers in order to store the semantic and structural information associated with a tetrahedral mesh of the myocardium : definition of endocardium, epicardium surfaces, scar regions and fibre orientation at each tetrahedron (see Figure 5). For the scene visualisation, SOFA integrated graphical interface based on 3D OpenGL was used. Volumetric meshes can be displayed with one dimensional textures either to represent the pseudo action potential or the isochrones. The multi-front anisotropic FMM was implemented as a SOFA solver component combined with the time-stepping mechanism already implemented in SOFA.

Benchmarks

We benchmarked the performance of our multifront anisotropic FMM (on a Intel Xeon CPU, 8 cores, 2.66 Ghz, 16 Go RAM) to evaluate the effect of mesh resolution and conduction anisotropy on our algorithm. In a first stage, we used a simple geometry representing a beam of size 5×2×2 cm with different mesh resolutions 6: Those tetrahedral meshes have been tessellated from a regular voxel grid. In those experiments, we assumed that the conduction velocity along the beam longest axis was 50 cm/s (typical velocity in cardiac tissue) starting from a single point located at a corner of the beam. In case of anisotropic simulation, the conduction velocity was set to be three times slower along the beam shortest axis. In Figure 1, isochrones of two propagations are displayed, showing the effect of anisotropy. The time step for those benchmarks has been set to 3 ms thus requiring around 33 to 43 iterations before depolarising the whole beam. In Figure 7 we plotted the computation time spent during each time step for the isotropic and anisotropic cases and for the three mesh resolution. As expected, in the isotropic case, the beam is depolarised in a bit more than 100 ms and the computation time reaches a plateau when the front corresponds to a section of the beam and therefore the same number of nodes needs to be updated at each time step. In the anisotropic case, such plateau does not exist since propagation fronts are no longer planar, having curved shapes varying along the beam. The wave required 130 ms to depolarise the whole beam since the wave has a slower transversal speed. Table 3 compares the total computation time (i.e. the area below each curve of Figure 7) as a function of mesh size and anisotropy. Clearly the total computation time is proportional to the number of mesh nodes and the anisotropic propagation requires in average 27 % more computation. It should also be noticed that the peak computation time per time step is also 27 % greater in the anisotropic case than in the isotropic case. This means that the increased computation time due to anisotropy is not only caused by a longer propagation but also by more computations per time step. A second set of benchmarks was performed on a myocardium mesh consisting of 10 213 vertices and 47 382 tetrahedra (see Section 4.1). In this case, anisotropic propagation was simulated with the same anisotropy ratio and conduction velocity as in the beam experiment. Figure 8 (Top) shows the plot of the computation time per time step as done previously on the beams as well as the percentage of nodes labelled as accepted (or depolarised) or refractory (in their refractory period). One can see that the computation time curve is bell-shaped which matches the pattern observed on anisotropic propagation for the beam. Furthermore, one can see that the proportion of depolarised nodes increases to 100 % when the myocardium is fully depolarised. Then, this proportion slowly decreases as more and more nodes are labelled as refractory and then at rest. One should notice that the change of states (from accepted to refractory to rest) requires a negligible computation time. Finally, in Figure 8 (Bottom) we investigate the discrepancy between the simulation time and the computation time. Ideally, one would expect that computation time for a single time step of 3 ms is exactly 3 ms. However, unlike classical PDE iterative solvers, our algorithm is asynchronous in the sense that the number of nodes processed at each time step may greatly varies. In this figure, we com-pare the cumulative computation time in red and orange with the ideal cumulative time in case of synchronicity (computation time is equal to simulation time) in blue. One can see that the computation time is always lagging behind the simulation time. This delay grows quickly during the depolarisation phase (where the anisotropic FMM is computed) but resides during the repolarisation and diastolic phase (where only the states are updated). In the case of isotropic propagation the total computed time for a single heart beat is almost equal to the total simulated time (near 0.9 s). However, the perceived propagation is slightly distorted since the depolarisation is artificially slowed down while the diastolic interval is accelerated. The anisotropic propagation is even more distorted with a 60 % increase in total computation times. This study shows that while our proposed algorithm is fairly close to real-time computation (up to a factor 1.9), it entails some significant alterations in the perceived propagation which must be corrected by a further speed-up.

Radio-Frequency Ablation Simulator

In addition to simulating cardiac electrophysiology, we propose in this paper to develop a simulator to train cardiologists in performing RF ablation. To be complete such simulators should also reproduce the catheter navigation performed under X-ray guidance which contributes to the difficulty in performing the intervention. Several research groups [START_REF] Dequidt | Interactive simulation of embolization coils: Modeling and experimental validation[END_REF][START_REF] Luboz | Real-time seldinger technique simulation in complex vascular models[END_REF] and commercial companies [33,34] are developing such endovascular navigation simulators, but we focus on the simulation of electrophysiology propagation and radio-frequency ablation which are suitable for realtime interaction.

Simulation Scene

As explained previously, multi-front anisotropic FMM is an efficient computational method to approximate the propagation of transmembrane potential in the heart. We can thus visualise that potential as a coloured texture. In order to illustrate the radio-frequency ablation simulation, we build a simulation scene corresponding to ischemic Ventricular Tachycardia (VT), an example of cardiac arrhythmia. In this pathology, the myocardium includes several necrotic regions or scars that appear after an infarct. Those scars can be imaged with a specific MR protocol, late enhancement MR, and are known to correspond to regions of very low conduction velocity. In the scene, we have decided to create synthetically a partially transmural isthmus between the two scars regions, which occurs frequently on patients suffering from VT [START_REF] Stevenson | Identification of reentry circuit sites during catheter mapping and radiofrequency ablation of ventricular tachycardia late after myocardial infarction[END_REF]. This isthmus, represented in red in Figure 9, is composed of peri-infarct tissues for which the propagation of electrical waves is slowed down. To prevent fast heart beats and re-entries, radiofrequency ablation has been shown to be an efficient therapy if performed optimally [START_REF] Hanno | Catheter ablation of multiple ventricular tachycardias after myocardial infarction guided by combined contact and noncontact mapping[END_REF]. During this procedure, catheters are used to permanently burn cardiac tissues with electromagnetic radiations having low voltage and frequencies in the radio-frequency range. When heated, the damaged tissue is no longer able to conduct electrical waves preventing the arrhythmia from occurring.

Concretely in SOFA, to simulate this pathology, a volumetric mesh of the right and left ventricles with scars has been generated from MRI. This mesh is composed of 10 213 vertices and 47 382 tetrahedra. Figure 9 shows the conductivity map on this mesh. Healthy tissues have a conduction velocity of 50 cm/s whereas the scar regions (in blue) have no conduction velocities and the isthmus (in red) has a conduction velocity of 30 cm/s.

Interactive Simulation

We first developed two mechanisms with which the user can interact with the scene. The first one consists in locally modifying the conductivity of a region of the myocardium whereas the second one performs an electrical stimulation which generates an additional electrophysiological wave by forcing the state of neighbouring nodes.

We bound those two interactions to mouse events and then, to improve the realism of the simulation, we interfaced them with a 3D user interface. Ideally, the user should interact with a force-feedback interventional radiology simulator which would allow the user to navigate a catheter through blood vessels and perform various actions (measuring potentials and ablating tissue). However, as a first prototype, we are using the Xitact ITP instrument [33] which is designed for laparoscopic simulation with force-feedback. This tool allows 3D positioning and consists of one tracked instrument with a button and two pedals.

For this instrument the position and orientation of its tips can be read through a serial port at a high frequency. A complete description of the handling of haptic device, positioning and force feedback in the SOFA framework, can be found in [START_REF] Saupin | Contact model for haptic medical simulation[END_REF].

Instrument Model

Due to the nature of the 3D user interface, we assume that catheters are rigid instruments thus avoiding the simulation of flexible catheters that would be used for this kind of intervention. Figure 10 shows that the instrument is modelled as a thin rigid shaft with a sphere at its extremity. This tip is used to simulate the catheter interaction (potential stimulation and ablation). The sphere serves as a collision model for the collision detection with triangles on the myocardium surface (see Figure 11). Thus, as soon as the sphere intersects a triangle, a collision is detected and the information is sent through the default collision pipeline of the SOFA framework [START_REF] Allard | Sofa an open source framework for medical simulation[END_REF]. Upon the detection of a collision, no force-feedback is sent back to the device in this first prototype. Instead, to guide the user, the part of the mesh being in collision is highlighted as a selected area and its scale may be modified as seen in Figure 11. Then, the user can push a button on the ITP device either to simulate cardiac tissue ablation (changing its conduction velocity to zero) or to perform an electrical stimulation. A foot pedal is used to change the nature of the interaction. In future work, we plan to use the Xitact VIST-C device for this simulation. By integrating our proposed algorithms of fast cardiac electrophysiological modelling with existing developments on catheter navigation already implemented in the SOFA framework [START_REF] Dequidt | Interactive simulation of embolization coils: Modeling and experimental validation[END_REF], we think that a more realistic simulator should be built in the near future, although some work remains to be done.

Simulation Results

To display the electrophysiology simulations, we use a colour map from blue to red, to represent the transmembrane potential of each vertex. Depolarised vertices are red, repolarised vertices are blue. At the initialisation stage, all vertices are set to the resting state.

Then, a stimulation is performed periodically at the apex of the right ventricle at a period of 300 ms to simulate a clinical procedure called VT Stim. This procedure is performed to test the inducibility of ventricular tachycardia of a patient and decide whether a defibrillator should be implanted. The electrical wave resulting from the stimulation takes approximately 250 ms to cover the two ventricles.

The first row of Figure 12 shows the resulting wave propagation (the cardiac fibres responsible for the anisotropy can be seen in the first image). As expected, the scars remain blue, since it is non conductive whereas the isthmus zone has a different behaviour than neighbouring healthy tissues. Indeed, the wave enters the isthmus through the upper side (A) and then exits by its lower bound (B). Then, due to different parameter values of the restitution curve, these vertices have a longer refractory period than healthy tissue (second row of Figure 12). As a result, when the second wave (wave T+1) reaches the upper bound of the isthmus, it can no longer go through this path. However, as the isthmus is no longer in refractory period when the second wave reaches the point B, it can then enter the isthmus through this point. This pattern is called an unidirectional block, and it can generate a re-entry wave (third row of Figure 12). In an ischemic VT, this re-entry wave can turn around the scars through the isthmus and then create additional ventricular contractions leading to tachycardia. In those cases, when the wave T+1 goes through the isthmus from the point B, on the other side (point A), healthy tissues are already in repolarised state. Thus, when reaching the point A, the wave will be able to propagate into the heart, as it is the case in the third row of Figure 12, and potentially sustain itself by repeating this pattern.

Additionally to the simulation described above, we simulated interactively the radio-frequency ablation therapy. In this case, the aim is to close the isthmus in order to prevent re-entry loops. With this set-ting, several radio-frequency ablations are performed around the isthmus in order to burn its tissue and get a large region with zero conduction velocity which blocks electrical waves. To ensure that the path is completely blocked, stimulation can be performed during the intervention by actioning a foot pedal. By applying stimulations near the scars, one creates a new potential wave to check that the isthmus has been successfully ablated. There is an associated video.

Discussion

To assess the accuracy of the simulation, it is possible to compare the simulated isochrone map with the one measured by catheters during the intervention using contact (i.e. the NavX or Carto systems) or non-contact (such as the Ensite balloon) endocardial mappings. Further evaluation can be performed by comparing stimuli conditions which can induce VT.

In terms of haptic simulation, no validation from clinicians can be performed at the moment since our prototype does not reproduce virtual catheter guidance to the heart, as it uses laparoscopic instruments for force feedback.

Conclusion and Perspectives

We presented a new algorithm to achieve close to real-time simulations of cardiac electrophysiol-ogy. This model opens up possibilities for filtering (smoothing) and interpolating electrophysiology recordings, which could provide diagnosis and therapy planning tools. An excellent example application is the planning of radio-frequency ablation. Such a real-time model can be embedded in an intervention simulator to test several ablation strategies.

To improve the model, future work includes speeding-up the algorithm to reach synchronicity and taking into account the impact of the curvature of the front on the propagation speed in the multi-front case. Accuracy of the method with respect to the mesh and the gradient approximation must also be carefully estimated.

This model has been integrated in the open source framework SOFA. Using a 3D user interface, we have presented a complete pipeline for the interactive simulation of cardiac radio-frequency ablation. A first prototype has been demonstrated for ischemic ventricular tachycardia where re-entries may occur through an isthmus surrounded by scar regions. A 3D rigid instrument allows the user to interactively simulate radio-frequency ablation by zeroing the conduction velocity around intersected tetrahedra. The user can also apply focal stimulations to check the effect of ablations.

There are several ways to improve the realism of the current prototype. For instance, the electrophysiology model could be extended to better represent the pathological behaviour and to include patientspecific tissue properties, such as conductivity and restitution curves. Real-time coupling between mechanical deformations and electrophysiology propagations should be developed in order to provide additional realism. Finally, catheter navigation should be made more realistic by using interventional radiology force-feedback simulators.

Annex

In order to compute the depolarisation time of a vertex X i within a tetrahedron (X i , X, Y, Z) , Algorithm 4 is presented. It invoves three alternative minimisations depending on the number of known neighbours among the three vertices X, Y , Z. These min-imisations can be computed analytically and correspond to three distinct cases [START_REF] Qian | A fast sweeping method for static convex hamilton-jacobi equations[END_REF]:

1. If X i has only one known neighbour X then the update does not include a minimisation and it is computed directly by linear extrapolation. 2. If X i has two known neighbours X and Y , then the minimisation corresponds to solving the following quadratic equation involving v(p), the point on the segment [X, Y ] located along the characteristic direction, p being its barycentric coordinate (p ∈ [0, 1]) :

T (X) -T (Y ) + v(p) ′ D -1 ( --→ XX i - --→ Y X i ) F [v(p)D -1 v(p)] = 0, with v(p) = --→ XX i p + --→ Y X i (1 -p)
3. If X i has three known neighbours X, Y , Z, the solution becomes more complicated and we do not include it here in the interest of space. We refer the interested reader to [START_REF] Qian | A fast sweeping method for static convex hamilton-jacobi equations[END_REF].

Algorithm 4 Computation of T (X i ) = U P DAT E(X i , X) in a Tetrahedrisation T (X i ) ← +∞ Update X i , the neighbours of X (which is known) within the surrounding tetrahedra T ET X Xi for all T ET (XX i Y Z) ∈ T ET X Xi = {T ET (XX i Y Z)|Y, Z ∈ N (X i )} do if Y, Z ∈ KN OW N then T (X i ) ← min{T (X i ), min (p,q)∈[0,1] 2 ([T (X)p + T (Y )(1p)]q + T (Z)(1q) + [v(p, q) T D -1 v(p, q)] 1/2 /F )} where v(p, q) = [ --→

XX i p + --→ Y X i (1 -p)]q + --→ ZX i (1 -q) else if Y ∈ KN OW N then T (X i ) ← min{T (X i ), min p∈[0,1] (T (X)p + T (Y )(1 -p) + [v(p) t D -1 v(p)] 1/2 /F )} where v(p) = --→ XX i p + --→ Y X i (1 -p) else T (X i ) ← min{T (X i ), T (X) + [ --→ XX i t D -1 --→ XX i ] 1/2 /F } end if end for

Figure 1 :

 1 Figure 1: Isochrones of the simulated front, red colour corresponds to small times, green to large times. (Top) Isotropic Fast Marching (Bottom) Anisotropic Fast Marching on tetrahedral unstructured grid, the anisotropic direction is vertical. Initialisation is in the front top corner.

Figure 2 :

 2 Figure 2: (Top) Impulse train with the Multi-Front Fast Marching (anisotropic case). Excited vertices are in red and refractory vertices are in green. (Bottom) Isosurface of the completely computed depolarisation front: a smooth front can be displayed instead of the pseudo potential.

Figure 3 :

 3 Figure 3: Transmembrane potential (red) computed from the discrete states (green).

Figure 4 :

 4 Figure 4: Restitution curves for isthmus (green) and healthy (red) tissues.

Figure 5 :

 5 Figure 5: Visualisation of the fibre directions stored in each tetrahedron of a volumetric mesh of the myocardium.

Figure 6 :

 6 Figure 6: Beam tetrahedral meshes having three different mesh resolutions.

Figure 7 :

 7 Figure 7: Computation time per time step as a function of mesh resolution and anisotropy.

Figure 8 :

 8 Figure 8: (Top) computational time per time step (in red) and percentage of point depolarised and in a refractory stage; (Bottom) Cumulative computation time for anisotropic and isotropic propagation for one heart beat compared to the cumulative simulation time.

Figure 9 :

 9 Figure 9: Conduction velocities of the myocardium. Zones in wireframe have normal conductivity while blue regions (scars) have lower values and the red region (isthmus) has intermediate values.

Figure 10 :

 10 Figure 10: Picture of the ITP Xitact tool and its corresponding model. The tool tip corresponds to a sphere collision model.

Figure 11 :

 11 Figure 11: Collision model representations. A sphere model is used for the tool whereas the myocardium collision model consists of tetrahedra. The contact region is shown in green.

Figure 12 :

 12 Figure 12: Simulation of the VT stim procedure for VT ablation. The left catheter is at the apex of the right ventricle to pace the heart at different frequencies. The right catheter is within the left ventricular cavity in order to ablate the isthmus between two scars responsible for re-entries. Electrophysiological wave propagation simulated with a multi-front eikonal model. Successive waves are represented as well as isthmus boundaries (T, T+1, A and B on first image). There is an associated video.

Figure 13 :

 13 Figure 13: Simulation of the ablation of the isthmus. Colours reveal the underlying conductivity (green: healthy, red: scar).There is an associated video.

Table 1 :

 1 Parameter values used for the restitution curve during the VT simulation.

Table 2 :

 2 Meshes used for the benchmarks

Table 3 :

 3

Total computation time to depolarise a beam as a function of mesh size and anisotropy ratio.


H Delingette

F Billet

K C L Wong

M Sermesant

K Rhode

M Ginks

C A Rinaldi

R Razavi

N Ayache

Personalization of Cardiac Motion and Contractility from Images using Variational Data Assimilation

Keywords:

Personalization is a key aspect of biophysical models in order to impact clinical practice. In this paper, we propose a personalization method of electromechanical models of the heart from cine MR images based on the adjoint method. After estimation of electrophysiological parameters, the cardiac motion is estimated based on a proactive electromechanical model. Then cardiac contractilities on two or three regions are estimated by minimizing the discrepancy between measured and simulation motion. Evaluation of the method on three patients with infarcted or dilated myocardium is provided.

I. INTRODUCTION

The objective of biophysical cardiac modeling is to improve the understanding of the cardio-vascular system by performing numerical simulation of the cardiac function from its mathematical description. Simulations of the whole organ have reached such a degree of realism that it is now possible to compare them quantitatively with available cardiac images and signals acquired routinely on patients. This has led to a new vision and a potential impact of those models in the clinical practice: improving the diagnosis of cardiac diseases and planning therapies (such as Cardiac Resynchronization Therapy [START_REF] Sermesant | Personalised electromechanical model of the heart for the prediction of the acute effects of cardiac resynchronisation therapy[END_REF]) from personalized biophysical models. Model personalization consists in optimizing some parameters of the model such that it behaves in adequacy with patient specific datasets (images and signals). This personalization can be split into two separate stages: geometrical personalization and biophysical personalization (see Fig. 1).

In the former case, the computational mesh is built from patient-specific anatomical images based on medical image analysis algorithms. In the latter case, the parameters of the models, initial and boundary conditions are optimized such that cardiac simulations match observations from images. Solving this inverse problem is a complex task for the following reasons: i) it is very time consuming and convergence is not guaranteed; ii) the information extracted from images is usually very sparse spatially and temporally iii) only a subset of the parameters may be observed and thus estimated. Furthermore, parameter estimation of dynamic systems often leads to mathematically involved methods.

Several approaches for cardiac model personalization have been proposed with various parameters to be optimized and optimization methods. In this article, we demonstrate a novel personalization method of a 3D cardiac model from catheterized electrophysiology data and cine MR images from three patients. The biophysical personalization is broken into three successive stages: electrophysiology, kinematics and mechanical personalization. A first contribution of our approach is to include an electromechanical model for the motion estimation stage (kinematics personalization) in order to produce a regularized motion. A second contribution is to propose an adjoint based method to estimate contractility parameters on the right ventricle (RV) and left ventricle (LV) as well as scar regions. Finally, we present a first automated method for biophysical personalization combining electrophysiology and mechanics.

II. ELECTROPHYSIOLOGY AND KINEMATICS PERSONALIZATION

A. Electrophysiology Personalization

In order to avoid accumulating sources of uncertainties, we are considering in this paper patient specific datasets that include a rich description of cardiac electrophysiology. More precisely, in addition to the acquisition of anatomical and cine-MR images, non-contact endocardial mappings have been acquired at the St Thomas' Hospital during electrophysiology studies in an XMR suite [START_REF] Rhode | A system for real-time XMR guided cardiovascular intervention[END_REF] both in sinus rhythm and during specific stimulation protocols. The extracted depolarization and repolarization isochrones then serve as input information to an electrophysiology personalization method [START_REF] Relan | Coupled personalisation of cardiac electrophysiology models for prediction of ischemic ventricular tachycardia[END_REF] which minimizes the discrepancy between measured and simulated isochrones. Its output is a set of global parameters and local parameters (electrical conductivities...) of the Mitchell-Schaeffer electrophysiology model which allow to interpolate, extrapolate and regularize the acquired isochrones. The personalized depolarization and repolarization times are then used to control the active contraction force of the mechanical model described in the next section.

B. Kinematics Personalization

Kinematics personalization consists in estimating the motion of cardiac structures from images. We use the same electromechanical model both to regularize the cardiac motion from cine-MR images and to estimate biophysical parameters. Indeed, this leads to an estimation of velocity fields which is consistent in terms of smoothness and spatial resolution with the ones produced by our electromechanical model used for the mechanical parameter estimation.

The kinematics personalization approach described in [START_REF] Billet | Cardiac motion recovery and boundary conditions estimation by coupling an electromechanical model and cine-MRI data[END_REF] is based on a proactive deformable model [START_REF] Sermesant | An electromechanical model of the heart for image analysis and simulation[END_REF] including three mechanical components: transmembrane potential propagation, active contraction forces, and passive biomechanics. The two ventricles are meshed with tetrahedra from anatomical MR images. The evolution of the displacement U of each mesh node is governed by the following equation:

M Ü + C U + KU = F b + F c + βF img (1)
where M, C, and K are the mass, damping, and stiffness matrices respectively, F b is the external load from boundary conditions (blood pressure, valves), F c is the force vector for active contraction. This contraction force is applied along a local fiber orientation and is controlled by several parameters including the maximum contraction parameter σ 0 . Furthermore, spring boundary conditions [START_REF] Sermesant | An electromechanical model of the heart for image analysis and simulation[END_REF] are applied to points located at the base of the ventricles.

In the kinematics personalization, we also add the term F img corresponding to a force vector which tracks salient image features in the image sequence. More precisely, F img is computed using a 3D block-matching algorithm to attract points towards the closest edge voxels. To balance between motion regularization and data attachment, we set the β parameters based on a sensitivity analysis [START_REF] Wong | Cardiac motion estimation using a proactive deformable model: evaluation and sensitivity analysis[END_REF]. Image forces are not physiology based since their sole purpose is to help tracking the cardiac motion. They are discarded during the mechanical personalization.

III. MECHANICAL PERSONALIZATION BASED ON VARIATIONAL DATA ASSIMILATION A. Adjoint Method for Parameter Estimation

Parameter estimation of biophysical cardiac models is a fairly new problem receiving a growing interest. Sequential and variational methods are among the most popular data assimilation techniques. Sequential methods such as rUKF have been proposed recently for the parameter estimation of contractility or passive stiffness parameters on synthetic data [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF], [START_REF] Xi | Myocardial transversely isotropic material parameter estimation from in-silico measurements based on reduced-order unscented kalman filter[END_REF]. While very promising, these methods do not handle well the phase changes of the cardiac cycle, which are an important aspect of the heart behavior. Also, preliminary work based on variational data assimilation (with [START_REF] Sundar | Biomechanically-constrained 4D estimation of myocardial motion[END_REF] or without [START_REF] Sermesant | Cardiac Function Estimation from MRI Using a Heart Model and Data Assimilation: Advances and Difficulties[END_REF] solving the adjoint equation) have also been performed on synthetic data.

In this article, we use a variational data assimilation approach based on the adjoint model and evaluate it on three cine MRI cases. It is based on the minimization of a functional measuring the discrepancy between simulation and measurements, in our case the node positions previously estimated from the kinematics personalization using image forces:

F = N n=0 i∈Vn 1 2 i Y n M -i Y n K 2
where N is the number of images in the image sequence, V n is the set of nodes visible in images n, i Y n M (resp. i Y n K ) is the position of the node i at time n of the mechanically (resp. kinematically) personalized mesh.

Among the various parameters of the model, we choose to estimate contractility parameters Σ 0 = {σ R 0 , σ L 0 , σ S 0 } of the right and left ventricles and eventually scar regions. Those parameters can be mostly observed during systole while passive stiffness parameters can be mostly observed during diastole. More precisely, to observe contractilities one would need to measure also the endocardial pressure in addition to cardiac motion. With only the latter, we can only estimate the ratio between contractilities and passive stiffnesses which can be considered as apparent contractilities.

The optimization of functional F(Σ 0 ) is performed using a quasi-Newton BFGS-B with bound constraint. This optimizer requires the computation of the functional gradient with respect to each contractility parameter σ 0 which is performed using the adjoint method. More precisely, the discretization of Equation 1 with finite elements (in space) and Houbolt integration scheme (in time) leads to a recurrence relation :

K t+1 Y t+1 M = A t 0 Y t M + A t 1 Y t-1 M + A t 2 Y t-2 M + F t where Y t M
is the node position vector at time t, K t and {A t i } are stiffness matrices and F t is the force vector at time t. Therefore minimizing functional F is equivalent to the minimization under constraint:

L({Y t M }, {P t }, Σ 0 ) = F({Y t M }, Σ 0 )+ t f t=0 (P t ) T K t+1 Y t+1 M - 2 i=0 A t i Y t-i M -F t
where P t is the Lagrange multiplier associated with the recurrence relation at time t and t f is the total number of iterations. In this functional, only the terms F t and eventually K t , A t depend on the contractility parameters Σ 0 . However, ∂L/∂Σ 0 also involves the Lagrange multiplier P t . By writing that ∂L/∂Y t M = 0, ∀t, one obtains a recurrence equation only involving P j , Y j and A j i j ∈ [t -3, t + 3] corresponding to the adjoint model.

In summary, to estimate ∂L/∂Σ 0 , we first compute the Lagrange multipliers using the recurrence equation found by setting ∂L/∂Y t = 0. This recurrence is solved backward starting from P t f = 0. We then compute ∂A t i /∂Σ 0 and ∂F t /∂Σ 0 and finally obtain ∂L/∂Σ 0 using chain rule derivation.

In our electromechanical model, ∂A t i /∂Σ 0 = 0 and all stiffness matrices are constant. However, we have to deal with the isovolumetric contraction and relaxation phases which provide additional constraints. Also, we removed the Windkessel model describing the hydraulic resistance and capacitance in the pulmonary artery and aorta in order to decrease the implementation complexity. This unfortunately also limits the realism in the simulation of late ejection and early filling phases. Our approach differs from [START_REF] Sundar | Biomechanically-constrained 4D estimation of myocardial motion[END_REF], where the adjoint method is considered in the continuous setting and then discretized which may add inconsistent discretization errors. The total computation times for one gradient estimation (with forward and backward simulations) is around 45 min on a regular PC, thus leading to computation times of 4 to 6 hours for the mechanical personalization.

B. Calibration Phase

As any gradient-based optimization, our optimization approach is sensitive to initial values. Therefore, we perform a calibration phase of the contractilities based on the ejection fraction. More precisely, kinematics personalization with a standard value of contractility (typically σ 6 10 4 P a) is first performed. This leads to a reasonable approximation of the cardiac motion from which the ejection fraction of both ventricles can be estimated. Then several simulations of the electromechanical model are performed (without image forces) with varying values of contractilities in both ventricles until an ejection fraction close to the estimated ones is obtained. We then restart the kinematics personalization with this calibrated set of parameters.

IV. EVALUATION ON SYNTHETIC IMAGES

We evaluate the performance of the joint kinematics and mechanics personalization method on synthetic data. Instead of adding noise to a known simulation, we designed a more comprehensive approach by generating synthetic image sequences using known motion and mechanical parameters. In this case, we considered the electrophysiology and cine MRI dataset from a patient with heart failure. A scar region was detected from late-enhancement MRI and a tetrahedral mesh including the scar region was created (see Fig. 2 (Middle)). We also used patient-specific endocardial activation map measured with the Ensite balloon (St Jude Medical), which was extrapolated to the myocardial volume using an electrophysiological model. The electromechanical model is then run with a fixed set of contractilities in the LV, RV and scar region. We then use the obtained motion field to generate a synthetic MRI sequence from the end-diastolic image.

We apply the calibration stage for both right and left ventricles to obtain a proper initial value of the contractilities (see Fig. 3 (Left)), to find that 10 5 P a is a good initial value. After kinematics personalization we obtain a mesh which serves as a reference for the mechanics personalization. After 7 iterations, the optimization stops and converges toward contractility values very close from their reference values, with an overestimation of the lowest contractility. On this synthetic sequence, we achieve a good motion recovery, with an average distance error of 0.15 mm between the personalized model and the motion tracking model and an average distance error of 0.14 mm (see Fig. 2 (Right)) with respect to the ground truth (used to simulate images).

V. EVALUATION ON CLINICAL DATA

This method was evaluated on the same data as the one used for the synthetic case, but with the actual cine MR image sequence. We estimated the apparent contractility parameters σ 0 for each of the three regions using the kinematics and mechanics personalization method. The calibration phase leads to an initial value of 10 5 P a for both ventricles. A good convergence of the 3 parameters is reached with a lower contractility parameter in the scar region as expected (see Fig. 4). The RV contractility is greater than the LV one which may be explained by a greater ejection fraction measured in the RV than in the LV. The optimization decreases the mean distance error with respect to tracked motion (distance between i Y n S and i Y n E ) over the whole cardiac cycle from 2.1 mm to 1.6 mm (see Fig. 5). The average error is comparable to the in-plane image voxel size (1.56 mm 2 ) but the maximum distance error was 9.7 mm which is significant. The region of largest errors is the base which may imply that the boundary conditions at the base should be improved. However, this is also the region with the largest longitudinal motion and thus with high uncertainty in motion tracking also caused by the large slice thickness (10 mm in such clinical routine sequences). The end diastole is the time at which the average distance error is maximum (2.32 mm). A quantitative comparison of motion with tagged MRI in short axis slices for this patient appears to indicate an excellent agreement for the radial motion (maximum errors of 1.2mm) but a decreased accuracy (average errors of 2.5mm) for circumferential motion probably due to the aperture problem. In order to evaluate the realism of the personalized model, we also compared the simulated pressure with the personalized model to the pressure which was invasively measured on that patient during the study. We obtained a fairly good match of this pressure during the contraction phase (see Fig. 4). The maximum value of the rise of pressure (dP/dt) max was decreased from 148 kPa/s to 103 kPa/s after personalization to be compared with a measured (dP/dt) max of 78 kPa/s. The personalization approach was also applied to two additional datasets (using also electrical mappings and cine-MRI) of patients suffering from dilated cardiomyopathy. The estimated contractilities for the RV and LV are rather low (see Fig. 5 (Right)), especially for patient 3. The adjoint method leads to a decrease of the average distance error to 2.33 mm for patient 2 and 2.19 mm for patient 3 (see Fig. 6) for images with respective isotropic voxel size of 1.44 mm and 1.52 mm. One limitation of this work is that it is restricted to the estimation of a limited number of contractility parameters defined on prescribed regions. One could iteratively refine the number of regions where parameters are estimated in a coarse to fine manner like in [START_REF] Relan | Coupled personalisation of cardiac electrophysiology models for prediction of ischemic ventricular tachycardia[END_REF].

VI. CONCLUSION

We presented a new method for the automated mechanical personalization of cardiac models and applied it to several clinical cases including electrophysiology and cine MRI data. Those results are encouraging and their analysis can help improving the model and the parameter set to be optimized. Those personalized models may serve to plan various therapies (e.g. Cardiac Resynchronization Therapy) by predicting the cardiac function after testing several therapeutic strategies in silico.

Fig. 1 .

 1 Fig. 1. Electromechanical Models Personalization from Images and Signals.

Fig. 2 .

 2 Fig. 2. (Left) Mesh and fiber orientation used for cardiac simulation; (Middle) Scar regions; (Right) Cardiac meshes on long axis image slice; green: with initial parameters; blue : after mechanics personalization; red : ground truth.

Fig. 3 .

 3 Fig. 3. (Left) Calibration curve of contractility for the left ventricle; (Right) Estimation of 3 regional contractilities.

Fig. 4 .

 4 Fig. 4. (Left) Estimation of the contractility parameters for 3 regions of the myocardium; (Right) Estimated and measured pressure curve.

Fig. 5 .

 5 Fig. 5. (Left) Cardiac meshes overlaid on image slice; green: with initial contractility parameters; blue : after mechanics personalization; red : target mesh (after kinematics personalization); movies are available as supplementary material; (Right) Contractility estimation for patients 2 and 3.

Fig. 6 .

 6 Fig. 6. Overlay of the cardiac mesh after kinematics personalization for patient 2 (Left) and patient 3 (Right) with the same convention as Fig 5; movies are available as supplementary material.