M Sermesant
email: maxime.serme-sant@inria.fr

H Delingette

N Ayache

I. INTRODUCTION

I N this paper, we introduce a new integrated three-dimen- sional (3-D) model of the left and right ventricles of the heart which can be used for the simulation and the analysis of cardiac pathologies. The overall principle is described in Fig. 1. Our in silico model includes knowledge coming from various disciplines including anatomy, electrophysiology and biomechanics, in a framework where it can be directly compared to in vivo measurements. By coupling this model to clinical data, one could simulate a number of pathologies or the effect of therapeutic actions, and extract a number of indexes of the cardiac function.

The computational modeling of the human body has been of increasing interest in the last decades [START_REF] Ayache | Computational Models for the Human Body, ser. Handbook of Numerical Analysis[END_REF], as it has benefited from progresses in biology, physics and computer science. It is now possible to combine in vivo observations, in vitro experiments and in silico simulations. There is an important literature on the functional imaging and modeling of the heart [START_REF] Katila | Functional Imaging and Modeling of the Heart (FIMH'01)[END_REF], [START_REF] Magnin | Functional Imaging and Modeling of the Heart (FIMH'03)[END_REF]. The following references provide examples of the measurement of electrical activity, deformation, flows, fiber orientation [START_REF] Macleod | Direct and inverse methods for cardiac mapping using multielectrode catheter measurements[END_REF]- [START_REF] Rhode | A system for real-time XMR guided cardiovascular intervention[END_REF], and of the modeling of the electrical and mechanical activity of the heart [START_REF] Mcculloch | Computational biology of the heart: from structure to function[END_REF]- [START_REF] Noble | Modeling the heart[END_REF]. Many of the functional models of the heart are direct computational models, designed to reproduce in a realistic manner the cardiac activity, often requiring high computational costs and the manual tuning of a very large set of parameters. In our approach, we rather select a level of modeling compatible with reasonable computing times and involving a limited number of parameters, thus allowing the potential future identification of the model parameters from clinical measurements on a specific patient (by solving the corresponding inverse problem). The first work in this direction was recently presented [START_REF] Sermesant | Progress toward model-based estimation of the cardiac electromechanical activity from ECG signals and 4D images[END_REF], [START_REF] Moreau-Villéger | Building maps of local apparent conductivity of the epicardium with a 2-D electrophysiological model of the heart[END_REF].

Also we would like this model to help the interpretation of cardiac image sequences. Cardiac image segmentation is still an active research area as reported in the survey by Frangi et al. [START_REF] Frangi | Three-dimensional modeling for functional analysis of cardiac images: a review[END_REF] and a special issue of IEEE TRANSACTIONS ON MEDICAL IMAGING [START_REF] Frangi | New Trends in Three-Dimensional Cardiac Image Analysis[END_REF]. The use of deformable models [START_REF] Mcinerney | Deformable models in medical images analysis: a survey[END_REF] is mainly limited to deformable surfaces [START_REF] Montagnat | A review of deformable surfaces: topology, geometry and deformation[END_REF], with an extension to spatial and time constraints [START_REF]4D deformable models with temporal constraints: application to 4d cardiac image segmentation[END_REF]. Whereas it ensures a better robustness against noise and can include trajectory constraints, there is no a priori knowledge introduced on the motion to help the segmentation. This was made possible with a four-dimensional statistical heart motion model computed from series of tagged MR images in [START_REF] Gérard | Efficient model-based quantification of left ventricle function in 3-D echography[END_REF]. This motion model allows a better initialization in the different images from a segmentation of the first image, thus, a better segmentation of the sequence. Another extended approach presented in [START_REF] Kaus | Automated segmentation of the left ventricle in cardiac MRI[END_REF] combines a point distribution model (PDM) and two coupled triangulated surfaces to segment the left ventricle, using also a motion model for the initialization. Deformable surfaces can also be coupled with prior segmentation, for instance multiscale fuzzy-clustering [START_REF] Sanchez-Ortiz | Automated 3-D echocardiography analysis compared with manual delineations and SPECT MUGA[END_REF].

Extensions of the PDMs used in segmentation methods are based on the active shape models (ASMs) and the active appearance models (AAM). A description and comparison of these two models can be found in [START_REF] Cootes | Statistical models of appearance for medical image analysis and computer vision[END_REF]. These models are now developed for spatio-temporal data. In [START_REF] Hamarneh | Deformable spatio-temporal shape models: extending ASM to 2D+time[END_REF], the ASM is extended to 2-D+time by introducing spatio-temporal shapes (ST-shapes). In [START_REF] Lelieveldt | Time continuous segmentation of cardiac image sequences using active appearance motion models[END_REF], the AAM framework is extended to 2-D+time by considering the image sequence as a single shape/intensity sample, giving the Active Appearance Motion Model (AAMM). These models can be theoretically extended to 3-D, but the size of the models and the difficulty to obtain good correspondences in 3-D images make it still a current research area.

At the same time, deformable templates evolved toward sophisticated approaches, for instance combined with Bayesian classification and Markov random fields [START_REF] Mignotte | Endocardial boundary estimation and tracking in echocardiographic images using deformable templates and markov random fields[END_REF], or coupling shape-space and Kalman-filter-based tracking [START_REF] Jacob | A shape-space-based approach to tracking mycordial borders and quantifying regional left-ventricular function applied in echocardiography[END_REF]. But few of these approaches integrate the volumetric aspect of human organs and the dynamic nature of the heart. Due to the complexity of the developed methods, it is mostly done by propagating the result obtained in one image (or one slice) to the next one [START_REF] Jolly | Segmentation of the left ventricle in cardiac MR images[END_REF].

Volumetric models have mostly been introduced in cardiac function analysis for interpretation [START_REF] Park | Analysis of left ventricular wall motion based on volumetric deformable models and MRI-SPAMM[END_REF], [START_REF] Park | A finite element model for functional aalysis of 4D cardiac-tagged MR images[END_REF] as they offer richer mechanical parameters. They were also introduced in deformation analysis for physically-based interpolation [START_REF] Ferrant | Registration of 3-D intraoperative MR images of the brain using a finite element biomechanical model[END_REF]- [START_REF] Azar | Methods for modeling and predicting mechanical deformations of the breast under external perturbations[END_REF]. We believe that volumetric models also allow one to introduce much more a priori knowledge on the organ directly in the segmentation process. It can be anatomical information, like fiber orientation, or mechanical behavior, to offer more reliable estimations of heart kinematics [START_REF] Papademetris | Estimation of 3D left ventricle deformation from echocardiography[END_REF]- [START_REF] Lo | H1 filtering and physical modeling for robust kinematics estimation[END_REF]. It can also be used to jointly estimate kinematics and mechanical properties of the myocardium [START_REF] Liu | Simultaneous estimation of left ventricular motion and material properties with maximum a posteriori strategy[END_REF]. Most of these approaches use hexahedral or tetrahedral meshes, but there are also alternative mesh-free methods proposed [START_REF]Meshfree representation and computation: applications to cardiac motion analysis[END_REF].

These models (geometrical and/or biomechanical) are passive models, i.e., they do not anticipate the cardiac motion, they only evolve under the action of 1) external image forces and 2) internal forces which constrain the regularity of the motion (geometrical models) or take into account the fiber orientations and a constitutive law (biomechanical models).

The key idea in this paper is to build a "ProActive Deformable Model" of the heart for image analysis. It is a volumetric deformable model of the heart integrating a priori knowledge on the motion in the segmentation process through the simulation of the electrical propagation and the mechanical contraction. In the classification proposed by Frangi et al. review [START_REF] Frangi | Three-dimensional modeling for functional analysis of cardiac images: a review[END_REF], the presented method would fit in the "continuous volumetric model" class. We believe that this new generation of physiology-based deformable models opens new possibilities in cardiac function analysis. Moreover, it allows us to introduce pathological prior information into image analysis, compared to statistical motion models built on volunteers described in the literature.

The proactive model we introduce here presents internal forces which create a complete contraction of the two ventricles synchronized with the electrocardiogram (ECG), therefore, the external forces only have to create local corrections to adjust the model to the boundaries observed in the cardiac images.

Using a model with physics-and physiology-based parameters one can simulate some cardiovascular pathologies and interventions. For instance, this could help devise techniques to make electrophysiology therapies shorter, less invasive and more successful.

The electromechanical model of the heart presented is based on mathematical systems of nonlinear partial differential equations, set on a 3-D domain, considering the ventricles as a continuum.

We present first the anatomical mesh construction, then the electrophysiology modeling and the contraction simulation, through an electromechanical coupling. The computed model is compared to measures from both the literature and medical images. Finally two applications are presented: pathology simulation and segmentation of a cardiac image sequences.

II. ANATOMICAL MODEL CONSTRUCTION

The myocardium is represented as a tetrahedral volumetric mesh including anatomical information. The process to build such a model is detailed in [START_REF] Sermesant | Deformable biomechanical models: application to 4D cardiac image analysis[END_REF]. The main anatomical information we use is the myocardium geometry, its division into different anatomical parts and the local orientation of the muscle fibers. It is difficult to obtain both realistic geometry and smooth fiber orientations in the same coordinate frame.

Thus, two models were built, coming from two different data sets on canine hearts. One comes from dissection data ("UCSD") measured in Auckland, New Zealand (P. Hunter group) [START_REF] Nielsen | Mathematical model of geometry and fibrous structure of the heart[END_REF] 1 with very smooth fiber orientations (see Fig. 4), due to the smoothing and the interpolation of the 256 original points done in the University of California, San Diego (A.

McCulloch group). 2 The other comes from diffusion tensor MRI ("DTI") acquired at Duke University (E. Hsu group) [START_REF] Hsu | Myocardial fiber orientation mapping using reduced encoding diffusion tensor imaging[END_REF], with a geometry closest to observed canine anatomies, for instance the right ventricle shape and the septum thickness (see Fig. 3), but noisier fiber orientations.

Both datasets have a resolution close to , which is small enough for our application, especially compared to the size of the mesh elements. Depending on the application, one or the other quality is preferred, thus guiding the model choice.

We acknowledge that the demonstration would be better with a whole human in vivo dataset, but diffusion tensor imaging is not yet possible in vivo, so we did our best to integrate the available data.

A. Volumetric Mesh Creation

The geometry can be extracted from different medical imaging modalities. From a 3-D image of the heart, the myocardium is segmented, using classical image processing methods like thresholding and mathematical morphology. Then, a triangulated surface of the myocardium is obtained using the marching cubes method [START_REF] Lorensen | Marching cubes: a high resolution 3D surface reconstruction algorithm[END_REF], and is decimated to the required size (typically 7 000 nodes for accurate simulation, or 1 500 nodes for the segmentation of cardiac images). Finally, a volumetric tetrahedral mesh is created from the triangulated shell, using the INRIA software GHS3D 3 (Fig. 2).

B. Anatomical Labeling

To better control and analyze the model during the simulation, we label the anatomical mesh into different regions. These regions were segmented in the myocardium from the Visible Human Project by Prof. Karl-Heinz Höhne group, Hamburg University [START_REF] Pommert | Creating a high-resolution spatial/symbolic model of the inner organs based on the visible human[END_REF]. This labeling is done by registering the mesh with the atlas image, and then assigning to each tetrahedron the main class corresponding to the voxels whose centers are lying inside this tetrahedron (these voxels are obtained by rasterization, the whole procedure is detailed in [START_REF] Sermesant | Deformable biomechanical models: application to 4D cardiac image analysis[END_REF]).

Fig. 3 presents the result of the anatomical regions assignment from the atlas to the mesh.

C. Myocardial Fiber Orientations

We present here the fiber orientations assigned to the two models from the two different datasets. The quantitative comparison of the behavior of the two models is beyond the scope of this paper, but the multiplication of available DTI data opens up possibilities to precisely analyze the influence and variability in cardiac fiber orientations. As diffusion tensor imaging is noisier near the surface of the myocardium, the fiber orientations are smoother in the wall than pictured here.

The knowledge of the myocardial fiber orientations (Fig. 4) plays an important role in the realistic modeling of the electrical and mechanical activity of the heart. Indeed, the conductivity is typically four times larger along the fibers than in the transverse direction, the orientation of the fibers creates a strong anisotropy in the constitutive law of the material, and also constrains the direction of the contraction stress.

III. ELECTROPHYSIOLOGY: TRANSMEMBRANE POTENTIAL SIMULATION

Many different models have been proposed to simulate the cardiac electrophysiology. They are divided into two main approaches.

• Biophysical or ionic models: cellular level simulation, using as variables the concentrations of the different types of ions, and integrating different ion channels based on Hodgkin-Huxley equations [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]- [START_REF] Noble | Improved guinea-pig ventricular cell model incorporating a diadic spaceI and I , and length and tension dependent processes[END_REF]. • Phenomenological models: more macroscopic models using a simpler system of equations to compute the cell potential without explicitly computing the concentrations of ion. It can be a bi-domain model, where the variables are the extra-cellular and intracellular potentials, or a mono-domain model where the variable is their difference, the transmembrane potential. One such model is the FitzHugh-Nagumo system [49]- [START_REF] Knudsen | Qualitative modeling of mechano-electrical feedback in a ventricular cell[END_REF]. As we model the electrophysiology mostly to control the contraction, we use the second approach, because the contraction is mainly related to the transmembrane potential. Moreover, for clinical use, only the extra-cellular potential can be measured, not the different ions concentrations, so we cannot adjust the parameters of the ionic approach from clinical data. 3 http://www-rocq.inria.fr/gamma/ghs3d/ghs.html

A. Transmembrane Potential

The transmembrane potential wave propagation is simulated using a system based on FitzHugh-Nagumo equations. This approach yields fast 3-D computations and allows us to capture the principal biological phenomena:

• a cell is activated only for a stimulus larger than a certain threshold; • the shape of the transmembrane potential does not depend on the stimulus; • there is a refractory period during which the cell cannot be excited; • any cell can be stimulated. Aliev and Panfilov developed a modified version of the FitzHugh-Nagumo equations adapted to the dynamics of the cardiac electrical potential [START_REF] Aliev | A simple two-variable model of cardiac excitation[END_REF] (1)

where is a normalized transmembrane potential and is a secondary variable for the repolarization. and control the repolarization, and the stimulation threshold and the reaction phenomenon. Throughout this manuscript, we use dots to represent partial derivatives with respect to time. The variable needs to be normalized in the Aliev and Panfilov equation in order to insure propagation (FitzHugh equation) and a proper coupling with the repolarization variable .

This model is simplified here: the term represents the influence of pacing frequency on the transmembrane potential duration and this property is not needed at the moment, so this term is neglected. Parameter values are derived from [START_REF] Aliev | A simple two-variable model of cardiac excitation[END_REF]: , , . To obtain the actual transmembrane potential in , we use the scaling . Similarly, time is normalized with an action potential duration (APD) of 1.0 in these equations. When used in the model, the time in the integration of this model is scaled to obtain a more realistic value . The orientation of the fibers is introduced through an anisotropic 3 3 conductivity tensor in an orthonormal basis whose first vector is along the local fiber orientation , with the conductivity in the fiber direction, and the conductivity anisotropy ratio in the transverse plane.

In Cartesian coordinates, it can be written:

, where denotes the fiber orientation, the tensor product (for a column vector ,

), and the identity matrix.

As previously mentioned, the conductivity in the fiber direction is typically four times larger than the conductivity in the transverse plane, therefore, a typical value of is . This yields a velocity of the propagation of the transmembrane potential typically two times faster in the fiber orientation than in the transverse plane (as the propagation speed of the transmembrane potential is proportional to the square root of the conductivity).

With an initial excitation above the threshold, the simulated transmembrane potential with this system is qualitatively similar to the transmembrane potential measured on cardiac cells (Fig. 6). For the sake of shape comparison, we present here the transmembrane potential of a frog cardiac muscle cell from [START_REF] Ruch | Physiology and Biophysics[END_REF] (the digital version is from [START_REF] Malmivuo | Bioelectromagnetism, Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF]).

B. Three-Dimensional Simulation of the Propagation

These equations are integrated on the 3-D volumetric mesh defined in Section II. A normalized transmembrane potential of 1.0 is imposed at the nodes corresponding to the Purkinje network terminations as an initial condition. The dynamic propagation can be represented by displaying an iso-surface of the transmembrane potential value. The complete propagation can be shown with the isochrones, where colors represent the different depolarization times (Fig. 5). The implementation of the model is described in Section V.

IV. BIOMECHANICS: ELECTROMECHANICAL CONSTITUTIVE LAW AND BOUNDARY CONDITIONS

A. Myocardium Mechanical Model

The myocardium is an active nonlinear anisotropic viscoelastic material. Its constitutive law is complex and must include an active element for contraction, controlled by the transmembrane potential computed in the previous section, and a passive element representing the mechanical elasticity. Several constitutive laws have been proposed in the literature [START_REF] Hunter | The analysis of cardiac function: a continuum approach[END_REF]- [START_REF] Caillerie | Toward a fiber-based constitutive law for the myocardium[END_REF]. These laws are designed to precisely fit rheological tests made on in vitro cardiac muscle.

Another approach is to model contraction from the nanomotors scale and build up a macroscopic constitutive law representing the phenomena encountered at the different scales, which is the approach followed by Bestel-Clément-Sorine [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF]. A detailed study of this complex model and one-dimensional (1-D) simulations can be found in [START_REF] Bestel | Modèle différentiel de la contraction musculaire contrôlée: Application au système cardio-vasculaire[END_REF] and [START_REF] Chapelle | A physiologically-based model for the active cardiac muscle contraction[END_REF]. This model is based on the Hill-Maxwell scheme, where muscles are represented by a combination of a contractile element, developing the stress tensor created by contraction, a series element, allowing isovolumetric contraction (especially in 1-D models), and a parallel element, mainly representing the elastic properties of the muscle.

B. Simplified Mechanical Model

The electromechanical model proposed here was motivated by the multiscale and phenomenological approach of [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF]. But it is specifically designed for cardiac image analysis and simulation. It is built in order to be computationally efficient and with few parameters, so we chose to simplify the constitutive law of [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF]. In our implementation, the model can be directly compared with in vivo measures through medical images. Despite its simplicity compared to other constitutive laws proposed in the literature, it reproduces quite well the global and local behavior of the myocardium.

The simplified mechanical model has the following components (see Fig. 7):

• an active contractile element which creates a stress tensor , controlled by the normalized transmembrane potential ; • a passive parallel element which is anisotropic linear viscoelastic and creates a stress tensor . The construction of these stress tensors from the transmembrane potential and the Lamé constants is detailed below.

For the electromechanical coupling, different laws have also been proposed [START_REF] Hunter | The analysis of cardiac function: a continuum approach[END_REF], [START_REF] Nash | Mechanics and material properties of the heart using an anatomically accurate mathematical model[END_REF]. We chose a simple ordinary differential equation to control the coupling, directly computing the contraction intensity from the normalized transmembrane potential. We believe that it is important to keep the model simple as relatively few clinical measures are available to adjust it. The contractile element is controlled by the normalized transmembrane potential through the ODE [START_REF] Katila | Functional Imaging and Modeling of the Heart (FIMH'01)[END_REF] with the time derivative of . As our normalized transmembrane potential is between 0 and 1 and the changes on depolarization and repolarization are abrupt, we can analytically approximate the solution of this equation by replacing with the value 0 or 1, using the current computed value thresholded at 0.5. It makes it possible to avoid time stepping the ODE and to directly control the parameters with the following coupling model as as [START_REF] Magnin | Functional Imaging and Modeling of the Heart (FIMH'03)[END_REF] with the depolarization time, the repolarization time, the heart period, the contraction rate, the relaxation rate, and . We added the constants and to better control the contraction stress increase and decrease. We can also add a time constant to and in (3) to model the delay between the electrical and the mechanical phenomena.

The 3-D contraction stress tensor is obtained with the formula , where denotes the fiber orientation and the tensor product. In the dynamics equation, when integrated over an element, it results in the force vector from Green-Ostrogradski formula, with the surface normal, and the element volume and surface, respectively. Contraction force is, thus, equivalent to a pressure applied along the fiber orientation.

This simplified constitutive law is represented by a damping matrix for the internal viscosity part, a stiffness matrix for the transverse anisotropic elastic part (parallel element) and a force vector computed from contraction (contractile element).

Once integrated into the dynamics equation, it writes (4) with the displacement vector, the diagonal mass matrix (mass lumping), a diagonal damping matrix, the anisotropic linear elastic stiffness matrix, the different external loads from the boundary conditions, and the force vector from the contraction.

As we consider the material linear elastic, but anisotropic, in small displacement formulation, is constant. The construction of the matrix is based on the finite element method with linear tetrahedral elements, with the derivation of displacements into the linearized strain tensor : and the Hookean constitutive law between Cauchy stress tensor and : , with the identity matrix, and , the Lamé constants. The details of implementation and the parameter values are in Section V.

The behavior of such a constitutive law is demonstrated on a cubic volume in Fig. 8. We can observe that the Lamé constants chosen to partly represent the incompressibility make the cube dilate vertically when it compresses horizontally.

C. Boundary Conditions: the Cardiac Phases

To simulate an entire cardiac cycle, the interaction of the myocardium with the blood is very important. This is why the different phases of the cardiac cycle have to be introduced, which implies different boundary conditions. The heart cycle can be divided in four phases: filling, isovolumetric contraction, ejection, and isovolumetric relaxation.

Four different boundary conditions are used on the mechanical model.

• Filling: a pressure is applied to the vertices of the endocardium. Its intensity is equal to the mean pressure of the atrium. It can be augmented during the P wave to introduce atrial contraction.

• Isovolumetric contraction: a penalty constraint is applied to the vertices of the endocardium to keep the ventricle volume constant. This penalty, which is iteratively estimated to counterbalance the contraction force, corresponds to the ventricular pressure during isovolumetric contraction. • Ejection: a pressure is applied to the vertices of the endocardium. Its intensity is equal to the mean pressure of the aorta (for the left ventricle) and the pulmonary artery (for the right ventricle). • Isovolumetric relaxation: a penalty constraint is applied to the vertices of the endocardium, keeping the volume constant. The penalty constraint is computed as follows: the volume at the beginning of an isovolumetric phase is computed, and at each iteration a pressure equal to (representing the ventricular pressure) is applied to the endocardial vertices, with the penalty factor and the current volume. Thus, if the volume is increasing, a negative pressure is applied, which tends to bring back the volume to its initial value. To ensure stability during this process despite the important stress developed, the time step has to be reduced during these phases (typically, it goes down from to ). In the current implementation, the atria pressures have two values (baseline and atrial contraction) and the arterial pressures (aortic and pulmonary) have a constant value.

To hold the mesh in space, we simulate the fibrous structure around the valves with springs having one extremity attached to a basal node and the other extremity attached to a fixed point.

To ensure mechanically smooth transitions between phases, change is automatically controlled in the following way.

• During filling, a pressure is applied to the endocardium. When the contraction starts, the contraction force will tend to eject blood, so when this force is more important than the applied pressure, the blood flow changes sign. As the blood is considered incompressible, the conservation of mass allows to compute blood flow directly with the ventricular volume time derivative. This is used to close the atrial-ventricular valves and start the isovolumetric contraction. • During the isovolumetric contraction, when the intensity of this penalty constraint is more important than the arterial pressure, the ventricular-arterial valves open, and the ejection phase starts. • During ejection, contraction force decreases after repolarization. When the flow changes sign, the ventricular-arterial valves close, starting the isovolumetric relaxation phase. • During isovolumetric relaxation, the penalty constraint represents the pressure, so when it is less important than atrial pressure, the atrial-ventricular valves open, starting the filling phase. Even with these completely independent conditions for the left and right parts of the heart, the two ventricles stay well synchronised, which shows that force development is coherent in the model. It allows us to adjust contractility parameters , and from the length of the different phases, and also from the atrial and arterial pressures (see Section VI-B1).

V. ELECTROMECHANICAL MODEL IMPLEMENTATION

The implementation of this model was done in C++, with a graphical interface in Tcl/Tk and OpenGL. It is ran on a Pentium PC, 2 GHz, and 1 Go of RAM. Parallel computations of the mechanical model are possible, which significantly decrease the execution time up to five processors, then the communication time becomes too important to achieve a real additional gain.

A. Electrophysiology Numerical Integration

The temporal integration is done with a fourth order Runge-Kutta scheme and the spatial integration is done with the Finite Element Method, using linear tetrahedral elements. The computation time step is and a 3-D simulation of the transmembrane potential during the cardiac cycle (0.85 s) takes around 5 mn on a standard PC with a 40 000 elements (7 000 nodes) tetrahedral mesh.

The parameters used in (1) are: , , , , and .

B. Biomechanics Numerical Integration

The mechanical model is integrated in time using the Houbolt semi-implicit scheme, and in space using the Finite Element Method with tetrahedral linear elements. Details of these methods can be found in many classical books, see [START_REF] Zienkiewicz | The Finite Element Method[END_REF] for instance.

We use the PETSc4 library for linear algebra operations, thus allowing distributed matrix storage, parallel preconditioning and parallel iterative solving. Details on the mesh partitioning, matrix assembly and parallel system solving are similar to [START_REF] Sermesant | Parallel implementation of volumetric biomechanical model and block-matching for fast nonrigid registration[END_REF].

To achieve this electromechanical simulation, we have to integrate two different phenomena: electrophysiology and biomechanics, and each of the models has a distinct inherent time step. If we call the electrical time step, the mechanical time step, and (respectively ) the current electrical (respectively mechanical) integrated time since the beginning of the simulation, at each instant of the cardiac cycle we integrate the least advanced phenomenon.

• If : we integrate the electrical phenomenon, and then , .

• If

: we integrate the mechanical phenomenon, and then , . The stability constraints from the boundary conditions are quite different during the different phases. We use an adaptive time step, with a time step times smaller during the isovolumetric phases.

The whole electromechanical cardiac cycle simulation with these boundary conditions takes less than 30 min on a standard PC (40 000 elements, 7 000 nodes, tetrahedral mesh). Half of the simulation time is devoted to the computation of the isovolumetric phases even if they represent only around 15% of the heart cycle, because they require a much smaller mechanical time step to achieve stability.

From the literature and the comparisons presented in next section, the following parameters are used:

(mass density), , , , , ,

. As no damping value could be found in the literature and its importance on cardiac function is still subject to debate, was used to control the numerical stability of the simulation and to obtain a reasonable visco-elastic behavior:

, with the identity matrix and the mass associated with the vertex. We use also mass lumping, the mass matrix is diagonal, and we associate the mass to each vertex corresponding to the accumulation of 1/4 of the mass of each tetrahedron containing .

VI. COMPARISON OF THE SIMULATED HEART CYCLE WITH MEASUREMENTS

A. Evaluation of the Transmembrane Potential Propagation

To simulate a realistic 3-D propagation of the transmembrane potential in the myocardium we need to determine the electrical onset for the initial conditions. The sinoatrial node is the natural pacemaker, located within the wall of the right atrium. It generates electrical impulses that are carried by special conducting tissue to the atrioventricular node. After reaching the atrioventricular node, located between the atria and ventricles, the electrical impulse goes down a conducting tissue (the bundle of His) that branches into pathways that supply the right and left ventricles. These paths are called the right bundle branch and left bundle branch respectively. The left bundle branch further divides into two subbranches (called fascicles). The extremities of these bundles are the Purkinje network, creating the junction between this special conducting system and the myocardium.

For our simulation, we need to locate these Purkinje network extremities, but it is hardly visible by dissection or by medical imaging. We used the measures from Durrer et al. [START_REF] Durrer | Total excitation of the isolated human heart[END_REF] which present illustrations of the isochrones in an isolated human heart, paced from the special conducting system (the version of these measures presented in Fig. 9 is from [START_REF] Malmivuo | Bioelectromagnetism, Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF]). The first isochrones in this article allow to visually locate the Purkinje network extremities on the endocardia of both left and right ventricles in the model and manually define them. Then the depolarization is simulated.

A first evaluation of the 3-D computation consists of comparing the resulting transmembrane potential isochrones with the measures from Durrer et al.. As we can see in Fig. 9, our simulation is qualitatively very close to the reported measures.

A more thorough evaluation of this electrophysiology model has been performed [START_REF] Moreau-Villéger | Building maps of local apparent conductivity of the epicardium with a 2-D electrophysiological model of the heart[END_REF] on canine hearts datasets from the Laboratory of Cardiac Energetics, National Heart, Lung, and Blood Institute, National Institutes of Health (NIH), showing that a local adjustment of conductivities could lead to a mean error in depolarization times of less than 5 ms (less than 5% error). Furthermore, it showed a good correlation between zones of low electrical conductivities and infarcted regions.

B. Evaluation of the Myocardium Contraction

We present in this section the comparison of local and global parameters of ventricles kinematics between our model and data extracted from medical images. The simulation results we present are stored after two simulated cycles, in order to obtain "natural" initial conditions from periodicity.

1) Volume of the Ventricles: We define a set of triangles representing the endocardium of each ventricle and we then close this surface with the barycentre of its edge to compute the inner volume of the ventricles (red surfaces and green lines in Fig. 10).

The evolution of the ventricle volume during the simulation of the cardiac cycle (Fig. 11) is very similar to the data available in the literature (see [START_REF] Silbernagl | Color Atlas of Physiology[END_REF] for instance). As we want to use this model for clinical applications, we also have to compare it with in vivo observations, i.e., medical imaging.

Automatic cardiac image segmentation is still a very challenging task, and manual segmentation of a full 4-D sequence is long and tedious. We present here a comparison with a volume curve extracted from a 3-D MRI sequence of a volunteer heart with the semi-automatic method detailed in [START_REF] Andriantsimiavona | Simple cardiac MRI segmentation[END_REF].

The ejection fraction computed from the simulated curve is 60%, compared to 63% from the measures. The main differences in the volume curve are during phase transitions, which are times when the intraventricular volume definition is not trivial in the images.

The evolution of these volumes makes it possible to adjust the contractility parameters: from the length of the isovolumetric contraction, from the length of the isovolumetric relaxation, and from the ejection fraction. The local motion described in the following sections results from these contractility parameters.

2) Local Apico-Basal Rotation: From the definition of the left ventricle endocardium, we can compute the inertia axis of the left ventricle (blue line in Fig. 10). We use this axis to compute the local apico-basal rotation of three points of the epicardium (Fig. 12) around this axis, throughout the cardiac cycle. This same rotation was measured for different points of the epicardium by Philips Research France through the analysis of tagged MRIs [START_REF] Allouche | A new kinetic modeling scheme for the human left ventricle wall motion with MR-tagging imaging[END_REF]. The values from the simulation show similar patterns and range to the measures (Fig. 13). Especially the opposite direction of rotation between the base and the apex is present both in the simulation and in the measures.

There are also some discrepancies between the two curves. For instance in the rising part of the apical rotation. This twisting motion originates from the fiber orientations but also from the isovolumetric phase and the activation sequence. It is still a phenomenon not completely understood, and as it is the result of many different elements, it is difficult to explain precisely these discrepancies. We will test the influence of different excitation sequences and fiber orientations on these curves to explore this.

The simulation of the transmembrane potential propagation and the inclusion of the different phases of the cardiac cycle is, thus, important to recover local parameters of the cardiac motion.

3) Local Radial Contraction: Another important local parameter of the cardiac function is the radial contraction, which measures the variation of the distance from a point to the central axis, throughout the cardiac cycle. The same inertia axis as for the rotation is used to compute the radial contraction of three points of the left ventricle endocardium (Fig. 14) during the simulated cycle. The same radial contraction was measured for different points of the myocardium by Philips Research France through the analysis of tagged MRIs [START_REF] Allouche | A new kinetic modeling scheme for the human left ventricle wall motion with MR-tagging imaging[END_REF]. The simulated radial contraction shows similar patterns with the measured one, in terms of range of values and profile (Fig. 15), which confirms the fact that the simulated ejection fraction is close to the real ones and that the model has a good local behavior. This radial contraction is responsible for the evolution of the wall thickness during the cardiac cycle, which is also a clinical index of the cardiac function.

VII. APPLICATION TO PATHOLOGY SIMULATION AND INTERVENTION PLANNING

Such a model enables the simulation of different cardiovascular pathologies, at the electrophysiological level or mechanical level. The observation of the consequences of these pathologies on the simulated cardiac function could help understand the phenomena, test different therapy strategies and plan interventions. A brief presentation on how this could be applied to different pathologies follows.

A. Ectopic Focus and Bundle Branch Block

An ectopic focus can be introduced by including an additional excitation point to the normal Purkinje extremities, with its own excitation sequence [Fig. 16 (left)]. A bundle branch block can be simulated by removing the Purkinje network extremities in one of the ventricles [Fig. 16 (right)].

B. Fibrillation

It has been shown that some cases of cardiac fibrillation are the result of a spiral of depolarization meandering in the myocardium. Such spirals can be simulated with the chosen model, using appropriate initial conditions like the wave-break method (Fig. 17) [START_REF] Krinsky | Spread of excitation in an inhomogeneous medium (state similar to cardiac fibrillation)[END_REF]. Studies of these spirals could help design more efficient defibrillators [START_REF] Gottwald | Spiral wave drift induced by stimulating wave trains[END_REF], by using a better defibrillation timing and, thus, less energy.

C. Radio-Frequency Ablation

The presented model was coupled with a force-feedback 3-D interface (Phantom from SensAble Technologies). It allows one to point locations on the 3-D model and change the local conductivity, whilst the simulation is running (Fig. 18). Although the simulation of the electrophysiology is not in real-time, the haptic device interface has to be, otherwise it would not be intuitive to control. We, thus, achieve an interactive change of the local conductivity: the ablated area is assigned a null conductivity in real-time, but the operator waits a few minutes for the completion of the simulation in order to observe the effect on the whole depolarization wave.

D. Infarcted Area

Some tissue pathologies, like infarcted areas, can be introduced in the potential propagation and in the mechanical contraction. Different effects can be investigated through simulations, for example the influence on the ejection fraction, which decreases from 65% to 55% in the simulated case (Fig. 19).

These are only generic simulations to present the idea of applying this modeling to pathology simulation. Validation patient data is in progress within the Cardiac MR Research Group, King's College London, Guy's Hospital, London [START_REF] Sermesant | Simulation of cardiac pathologies using an electromechanical biventricular model and XMR interventional imaging[END_REF]. 5 VIII. APPLICATION TO IMAGE SEGMENTATION: A PROACTIVE DEFORMABLE MODEL One of the applications of this model is for cardiac image segmentation. The key idea is to build a "Pro-Active Deformable Model" of the heart for image analysis, i.e., a volumetric deformable model of the heart integrating a priori knowledge on the motion, through the simulation of the electromechanical contraction. The internal forces regularising the deformation are computed from the electromechanical model previously presented and the external forces are computed from the image features. We, thus, solve the new dynamics equation introducing these forces [START_REF] Faris | Novel technique for cardiac electromechanical mapping with magnetic resonance imaging tagging and an epicardial electrode sock[END_REF] with the contraction forces, the image forces, and the weighting parameter for the image forces.

A. Internal Forces

The internal forces are computed from the electromechanical rheological model, thus, introducing the simulated contraction. We use the model presented in Section III-A to compute the action potential propagation and the model of Section IV-B to compute the mechanical contraction. For the time synchronization, information on the image sequence acquisition allows us to know the heart beat duration, the R wave position and the timing of each image in the cycle. This is used to trigger the transmembrane potential propagation and adjust the transmembrane potential duration, as well as to compute the external forces.

For segmentation purposes, some boundary conditions, like pressure and isovolumetric phases, are partly included in the image information. The added stiffness to represent the valves is also a part of the image information. This is the reason why no mechanical boundary conditions other than image forces are applied when the model is used in this image segmentation framework.

Usually, segmentation methods use one reference position of the model per image, corresponding to the previous image final position. We believe it is important not to reset the strain and stress at each time frame in order to capture all the characteristics of the deformation of the myocardium. Thus, we want to use only one reference position of the heart for the complete sequence.

Due to the large difference in shape between the end-diastolic position and the end-systolic position, using a priori knowledge on the motion through the contraction simulation helps to recover this deformation, as the image forces only have to correct the predicted deformation, not create it from the end-diastolic shape.

B. External Forces

The external forces are introduced as a load applied to the mechanical model. For each surface node of the mesh, we look for a corresponding boundary point in the image voxels lying along the surface normal, a classical approach for deformable models in computer vision, see Fig. 20. The boundary point is selected among these voxels from intensity, gradient direction and gradient value criteria [START_REF] Montagnat | Anisotropic filtering for model-based segmentation of 4D cylindrical echocardiographic images[END_REF]. Then a force is applied to this node, proportional to the distance to this boundary voxel and oriented in its direction.

These forces can be different on each of the different anatomical regions of the model, depending on what is visible in the different parts of the image and on the image intensity characteristics in these regions. If some parts of the myocardium are not visible, the external forces can be removed for the corresponding regions and only the internal forces will make these vertices move.

C. Global Adjustment to Patient Anatomy

The correspondences between the surface nodes and the boundary voxels can be used to globally adjust the mesh to the patient anatomy. After a rough alignment of the ventricles, we compute iteratively the best rigid transformation between the mesh and the image. After convergence, we compute the best similarity and finally the best affine transform. We use a new criterion proposed by X. Pennec, which is symmetric and "invariant" with respect to the action of an affine transformation to both the model and the image data. This feature helps to avoid convergence toward a singular affine transformation (null determinant) when the number of matchings is low, which often happens when fitting to noisy images. Let and be the matched model and image data points, the affine transformation and the translation. The criterion to minimize is with the identity matrix.

This criterion lends to an adequate initialization of the model with a global transformation even in noisy images. Details on the construction of this criterion and on the closed form solution can be found in [START_REF] Sermesant | Deformable biomechanical models: application to 4D cardiac image analysis[END_REF]. The model can then be fitted with a better accuracy using local deformations.

D. Sequence Segmentation: Time-continuous Image Force Field

As we use an electromechanical model for the internal forces, we need to integrate it according to time steps given by stability constraints, which are independent from the image acquisition time resolution. We create a "time-continuous image force field" by using the two images of the sequence surrounding the current integrated time in the cycle to interpolate the force to apply from the two forces computed within each of these images. This ensures a smooth evolution of the mechanical boundary conditions, thus, a better stability and segmentation process.

If we have an image at instant and the next one is at instant and the segmentation process is at time in the cardiac cycle, between and , then the applied image force is

E. Results

As a feasibility study of introducing contraction simulation in image analysis, we chose a modality where the definition of the external forces would be rather straightforward, in order to emphasise the influence of the proactive internal forces. Thus, we use SPECT imaging where boundary definition is quite clear to demonstrate the effect of the active internal forces, without too much influence of the choice of external forces. In other modalities, it is often difficult to obtain robust and consistent external forces computation, which is one of the motivations for introducing more prior knowledge, but also a major drawback when one want to compare different internal forces.

We present here the left ventricle volume curves obtained when segmenting a SPECT image sequence with a passive biomechanical deformable model and with the proactive deformable model, see Fig. 21. They are compared to reference values obtained with a deformable surface (simplex mesh) semi-automatically adjusted to each image of the sequence, without any temporal continuity. The balance between internal and external forces was optimized to obtain the best possible result with each model.

The consequence of the introduction of a priori knowledge on the motion in the segmentation is a better estimation of the volume, especially of the end-systolic position (see Figs. 22 and23). The ejection fraction, which is an important clinical index of the cardiac function, is then more accurately computed. The ejection fraction computed from the reference values is 66%. We obtain an ejection fraction of 68% with the proactive model (the electromechanical model alone has an ejection fraction of 62%) whereas the passive model results in only 53%.

Furthermore, the needed weight of the image forces, in (5), is ten times smaller with the proactive model, because it only has to correct the predicted motion, not to create the whole deformation from the end-diastolic position.

The model evolution gives a continuous estimation of the myocardium position throughout the cardiac cycle. It allows us to interpolates image information to correspond best at each time when an image is available and to continuously deform in between.

Moreover, using this proactive model gives a priori information on the local tangential motion (torsion) which is hardly visible in current medical images (without using tags). This kind of model could help recover this motion, which is also quite important in cardiac function.

The segmentation of a full image sequence with the proactive model takes less than 5 min on a standard PC. We use a coarser mesh (8 000 elements, 1 500 nodes) which is better adapted to the limited image resolution. The isovolumetric phases are not simulated in this case, saving a significant computational cost. The whole segmentation time does not depend on the time resolution of the image sequence as only the electromechanical phenomena control the integration time steps, and then image forces are computed for each of these integrated time steps.

IX. CONCLUSION AND PERSPECTIVES

The design of computational models of human organs is a new research field which opens new possibilities for medical image analysis and therapy simulation: this article presented a number of steps toward this goal in cardiac imagery, and must be understood as a preliminary proof of concept in this research direction.

The model we presented was designed at a macroscopic level with a limited number of internal parameters. Given the high complexity of cardiac motion, composed of different twisting rotations and radial and axial contractions, the proposed model still allowed realistic simulation of the heart motion while allowing reasonable computing time. We also showed that our "proactive" deformable model using its internal contraction forces synchronized on the ECG was able to better recover the segmentation of the heart ventricles from a time series of cardiac images than a more classical "passive" deformable model which would deform under the action of image forces only.

We acknowledge the limitations of the proposed model, at the three levels of representation: anatomy, electrophysiology and biomechanics. Possible improvements of the model would include the integration of more anatomical structures (valves, atria), a more realistic electrophysiology model (biophysical models) and a more complex constitutive law. However, our objective is not to build the more complex and faithful heart model ever. Instead, we want to adapt the complexity of our model to that of the available observations (images, signals) in order to establish new tools for diagnosis, therapy planning and therapy guidance. Thus, the described model is a starting point that could be refined in complexity if required by a given clinical application.

We already plan to test some refinements. First of all, the problems created by the different properties of the available anatomical data motivated us for the definition of an analytical model of the myocardium (from ellipsoids) which would help in the coherence of the data used. We hope that the rapid expansion of the possibilities of medical imaging will make it possible to have a precise idea of the variability of all the anatomical data used, in order to be able to estimate the validity of using such generic data.

We will look into introducing the atria and proximal portions of arteries, to help define boundary conditions and extend the number of pathologies that can be simulated. We also want to introduce a Windkessel model, frequently used to link blood pressure and blood flow in the arteries in order to compute the after-load [START_REF] Sagawa | Translation of Otto Frank's paper "Die grundform des arteriellen pulses[END_REF], [START_REF] Stergiopulos | Total arterial inertance as the fourth element of the windkessel model[END_REF].

We will also reconsider the simplifications made in the constitutive law, both in the passive and active parts, to be able to simulate more realistic local stress. We plan to test exponential laws of the passive myocardium and integrate an influence of the deformation in the active part, to better represent known behavior. However, we need to quantify the main source of error, for instance with a sensitivity analysis. Both the simplification of the constitutive law and the definition of the boundary conditions have a great impact on the simulated motion. For the second point, we are also looking into MR imaging of blood flows to help better define boundary conditions.

Finally, if needed for some pathologies, the electrophysiology model would be refined, for instance with a bi-domain or ionic model. But all these evolutions will be balanced against the computational cost and the parameter number increase to still be compatible with clinical applications.

To achieve a complete validation and allow a quantitative use of such models, one has to design parameter estimation procedures. It would enable the automatic identifications of the internal parameters from a number of clinical measurements, in particular new imaging technologies of cardiac electromechanical activity [START_REF] Rhode | A system for real-time XMR guided cardiovascular intervention[END_REF]. Results on global adjustment of models from these measurements are very promising [START_REF] Sermesant | Simulation of cardiac pathologies using an electromechanical biventricular model and XMR interventional imaging[END_REF]. It is of special interest to analyze sensitivity and observability of these parameters. Ultimately, we would replace the current image forces which are not physiology-based with an identification procedure which would optimize the physical parameters of the model to match the image data.

An

  Electromechanical Model of the Heart for Image Analysis and Simulation M. Sermesant* , H. Delingette, and N. Ayache Abstract-This paper presents a new three-dimensional electromechanical model of the two cardiac ventricles designed both for the simulation of their electrical and mechanical activity, and for the segmentation of time series of medical images. First, we present the volumetric biomechanical models built. Then the transmembrane potential propagation is simulated, based on FitzHugh-Nagumo reaction-diffusion equations. The myocardium contraction is modeled through a constitutive law including an electromechanical coupling. Simulation of a cardiac cycle, with boundary conditions representing blood pressure and volume constraints, leads to the correct estimation of global and local parameters of the cardiac function. This model enables the introduction of pathologies and the simulation of electrophysiology interventions. Moreover, it can be used for cardiac image analysis. A new proactive deformable model of the heart is introduced to segment the two ventricles in time series of cardiac images. Preliminary results indicate that this proactive model, which integrates a priori knowledge on the cardiac anatomy and on its dynamical behavior, can improve the accuracy and robustness of the extraction of functional parameters from cardiac images even in the presence of noisy or sparse data. Such a model also allows the simulation of cardiovascular pathologies in order to test therapy strategies and to plan interventions. Index Terms-Cardiac image analysis, cardiac modeling, deformable model, electromechanical coupling, simulation of cardiac pathologies.

Fig. 1 .

 1 Fig. 1. Overview: electromechanical model of the heart, interaction with patient clinical data and applications: cardiac function analysis and simulation of cardiac activity for pathology simulation.

Fig. 2 .

 2 Fig. 2. Tetrahedral mesh of the bi-ventricular myocardium (40 000 elements, 7 000 nodes) from the UCSD data.

Fig. 3 .

 3 Fig. 3. Anatomical regions obtained from the visible human atlas with the DTI geometry: basal left endocardial ventricle (A), basal septum (B), dorsobasal left epicardial ventricle (C), basal right ventricle (D), basal left epicardial ventricle (E), apical right ventricle (F), apical left epicardial ventricle (G).

Fig. 4 .

 4 Fig. 4. Fiber orientations assigned to the myocardium mesh and elevation angle from the data interpolated in UCSD (left) and from the DTI (right). Blue and red colors represent vertical fibers and green represents horizontal fibers.

Fig. 5 .

 5 Fig. 5. (left) Iso-surface of the simulated transmembrane potential value, representing the propagation front, at one instant of the cardiac cycle on the UCSD geometry. The red side is the depolarized one. (right) Resulting isochrones after complete myocardium depolarization.

Fig. 6 .

 6 Fig. 6. (Top) Transmembrane potential simulated with simplified Aliev and Panfilov model. Time is normalized to obtain an action potential duration equal to 1 time unit. (Bottom) Measured transmembrane potential on a frog cardiac muscle cell.

Fig. 7 .

 7 Fig. 7. Scheme of the simplified rheological model with a passive elastic element E and an active contractile element E .

Fig. 8 .

 8 Fig. 8. Contraction (first row) and relaxation (second row) simulation on a cube. Fibers are horizontal and an initial transmembrane potential is applied to the left facet of the cube. Color represents the transmembrane potential value (dark blue: polarized, light red: depolarized).

Fig. 9 .

 9 Fig. 9. Transmembrane potential isochrones (in ms) measured by Durrer et. al (top row) compared with the simulated ones (bottom row), on the UCSD geometry.

Fig. 10 .

 10 Fig. 10. Triangles sets used to define the endocardia of each ventricle, on the UCSD geometry. The barycentres of the edges of these sets are used to close the ventricles in the volume computation.

Fig. 11 .

 11 Fig. 11. Measured left ventricle volume from MRI compared to the simulated cycle. The simulated volume values and the ejection fraction (60%) are similar to the ones measured in volunteer data (63%).

Fig. 12 .

 12 Fig. 12. Basal, equatorial and apical points on the epicardium used to observe the local rotation during the simulated cardiac cycle, on the UCSD geometry.

Fig. 13 .

 13 Fig. 13. (left) Twisting angle during the simulated cycle, for different points of the epicardium: base (black), equator (dashed blue) and apex (red). (right) Twisting angle extracted from tagged MRI by Philips Research France for different points of the myocardium.

Fig. 14 .

 14 Fig. 14. Basal, equatorial and apical points of the left ventricle endocardium used to observe the local radial contraction during the simulated cardiac cycle, on the UCSD geometry.

Fig. 15 .

 15 Fig. 15. (left) Radial contraction during the simulated cycle for different points of the endocardium: base (black), equator (blue 0 0) and apex (red). (right) Radial contraction extracted from tagged MRI by Philips Research France for different points of the myocardium.

Fig. 16 .

 16 Fig. 16. Electrophysiology pathologies simulation, presented with resulting isochrones, on the UCSD geometry. (left) Ectopic focus (part of a Wolff-Parkinson-White syndrome simulation). (right) Right branch block simulation.

Fig. 17 .

 17 Fig. 17. Simulation of a reentry spiral using a wave-break. (Top) Description of the wave-break method. (Bottom) Simulation on a cube of myocardium model. Color represents the transmembrane potential (light red: depolarized, dark blue: re polarized).

Fig. 18 .

 18 Fig. 18. (left) Phantom 3-D interface, from SensAble Technologies. (right) Radio-frequency ablation simulation, by modifying the conducting parameters of the model where it has been in contact with the tool (homogeneous grey areas), on the DTI geometry.

Fig. 19 .

 19 Fig. 19. Basal left epicardial infarcted zone (blue, on the UCSD geometry). No conduction and no contraction in this area. (right) Corresponding volume curves without (dashed) and with (solid) infarct. Ejection fraction decreases from 65% to 55%.

  5 http://www-ipg.umds.ac.uk/m.sermesant/index.php

Fig. 20 .

 20 Fig. 20. External forces computation. (left) Scan-line algorithm to extract the image voxels from a vertex P along the normal ñ to the mesh surface. (right) Region and boundary criteria on the extracted voxels to determine the boundary point corresponding to the surface vertex P .

Fig. 21 .

 21 Fig. 21. Segmentation of cardiac image sequence with the proactive model of the heart. The model is displayed with slices of one image of the sequence, but the external forces are computed with two images. The transmembrane potential value is color-coded on the model surface and the red segments represent the external forces on the surface vertices toward the corresponding images voxels.

Fig. 22 .

 22 Fig. 22. Comparison between the passive (black-cross), electromechanical (dashed) and proactive (grey) models. The reference values are the black squares. The combined use of electromechanical model with image information achieves a better segmentation of the image sequence.

Fig. 23 .

 23 Fig. 23. Comparison of the segmentation of a SPECT image sequence, in three orthogonal slices of the end-systolic position (red: passive biomechanical model, yellow: proactive model). The electromechanical model stays closer to the image boundary, especially near the base (axial contraction) and we can observe the right ventricle contraction of the model, even if the right ventricle is not visible in the image (due to a previous region of interest extraction).

http://www.bioeng.auckland.ac.nz/home/home.php.

http://cmrg.ucsd.edu/

http://www-unix.mcs.anl.gov/petsc/petsc-2/

ACKNOWLEDGMENT

The authors would like to thank for their collaboration Dr. E. Hsu from Duke University; Prof. A. McCulloch from University of California, San Diego; Prof. K.-H. Höhne and his group from Hamburg University; Dr. O. Gérard and Philips Medical Systems Research Paris; the Cardiac MR Research Group in Guy's Hospital, London and all the co-workers of the ICEMA collaborative research actions (http://www-rocq.inria.fr/who/Frederique.Clement/icema.html; http://www-rocq.inria.fr/sosso/ icema2/icema2.html) funded by INRIA; in particular, M. Sorine, F. Clément, D. Chapelle, J. Sainte-Marie, Y. Coudière, S. Lantéri, and Z. Li. The authors would like to thank Timothy Carter, Olivier Clatz and Philippe Moireau for thoughtful comments, and Rado Andriantsimiavona for the left ventricle volume segmentation from cardiac MRI. Color version of the presented figures and video animations are available on the Internet at: http://www-sop.inria.fr/epidaure/personnel/Maxime.Sermesant/gallery.php