Abstract
In this work, we develop an interactive framework for rehearsal of and training in cardiac catheter ablation, and for planning cardiac resynchronization therapy. To this end, an interactive and realtime electrophysiology model of the heart is developed to fit patientspecific data. The proposed interactive framework relies on two main contributions. First, an efficient implementation of cardiac electrophysiology is proposed, using the latest graphics processing unit computing techniques. Second, a mechanical simulation is then coupled to the electrophysiological signals to produce realistic motion of the heart. We demonstrate that pathological mechanical and electrophysiological behaviour can be simulated.
1. Introduction
Cardiovascular diseases remain one of the most important causes of death worldwide. Cardiac arrhythmia and heart failure are lifethreatening cardiac diseases that consist of an abnormal electrical and/or mechanical activity of the heart muscle (myocardium). Such pathologies can occur upon changes in the heart structure following coronary artery disease (heart attack) or as chronic consequences of hypertension, diabetes or cardiomyopathy. Several therapies may be pursued depending on the pathology. For ventricular tachycardia, radiofrequency ablation of the ventricles may be performed, whereas cases of severe heart failure may be treated through cardiac resynchronization therapy (CRT).
Most of these procedures are minimally invasive and complex. They need to be performed by experienced and highly skilled cardiologists. Therefore, building training systems for rehearsing these endovascular interventions would help not only junior electrophysiologists to train for those procedures but also experienced electrophysiologists to rehearse complex procedures. There exists a limited number of commercial endovascular simulators [1,2], but they are somewhat limited because they rely on prestored electrophysiology data but not on electromechanical modelling (figure 1).
The objective of this paper is to propose an innovative framework for the interactive simulation of a cardiac electromechanical model. With the constant improvement in computational methods [3], patientspecific cardiac modelling is now being considered as a promising means of assessing therapy and increasing pathophysiological understanding. Physiological models not only can reproduce cardiac motion and electrophysiology signals but also can predict in silico the impact of therapies [4].
Our framework is mainly composed of two physiological models. The first is the modelling of cardiac electrophysiology. The electrical stimulation inside the myocardium is initiated by a natural pacemaker, called the sinoatrial node, located inside the right atrium. This generates a depolarization wave that propagates inside the heart wall, and generated by the exchange of ions through the membrane of cardiac cells. The difference in potential between the cardiac cells and the extracellular matrix is called the action potential (AP) or transmembrane potential. It drives the contraction of the cardiac fibres. Several mathematical models have been investigated by many research groups and they can be sorted into three different classes:

— biophysical models [5] are complex models simulating electrophysiology at the organ scale and involve a large number of parameters;

— phenomenological models [6–9] are simplified models (involving fewer parameters) derived from the biophysical models and capturing the AP shape at the organ scale;

— Eikonal models [10] are static nonlinear partial differential equations for the depolarization times obtained from the previous models; they are based on nonphysiological parameters and do not take into account all the complex physiological phenomena, such as reentries.
The second physiological model of our framework represents the mechanical coupling, i.e. how the depolarization of a cardiac cell generates tension in the sarcomeres, the basic unit of contraction, followed by relaxation. Such cellular mechanisms are highly complex and involve potentially many ionic types and parameters. Several models [11–15] have been proposed at the tissue scale to describe the relationship between the AP and active tension along the cardiac fibres.
In this study, we use the Mitchell–Schaeffer model [9] to describe the cardiac electrophysiology, because it has only two variables and depends on only a few parameters while capturing the main properties of restitution of the action potential duration (APD) and conduction velocities. We also chose the Bestel–Clément–Sorine (BCS) model to simulate electromechanical coupling owing to its compact nature and its consistency with the law of thermodynamics (balance of energy). The advantage in having a compact mechanical model, i.e. a model with few variables and parameters, lies in the ability to calibrate its parameters from a limited set of patientspecific data as shown by Marchesseau et al. [16] for cardiac mechanics. However, despite their relative simplicity, it must be stressed that both electrical and mechanical models lead to partial differential equations that are known to be timeconsuming to solve. For electrophysiological models, the steep potential upstroke constrains the time resolution, whereas, for mechanical models, isovolumetric constraints make the numerical system highly nonlinear.
Our main contribution lies in adapting and optimizing these models and the solver algorithms in order to reach a fully interactive electromechanical simulation of the heart. In the remainder, we present a realtime simulation of the ventricular electrophysiology that is compatible with user interaction. Then, an interactive mechanical model coupled with electrophysiology is computed to simulate the associated heart motion.
This study is therefore a first step towards realistic training systems and safer learning procedures for the future. It has been developed using SOFA^{1} [17], an open source framework for interactive numerical simulations in medicine.
2. Interactive cardiac electrophysiology
2.1. Geometrical and electrophysiology model
In this section, we first focus on the reconstruction of the biventricular geometry and on the simulation of electrophysiology based on the Mitchell–Schaeffer [9] model.
2.1.1. Geometry acquisition
We apply our interactive framework on both steadystate free precession (SSFP) and cineMRI sequences, which are, respectively, dedicated to anatomical description and motion tracking.
In order to build a patientspecific geometry from these sequences, a preliminary imageprocessing stage is necessary. First, we extract the two ventricles from the SSFP sequence, using a semiautomatic segmentation algorithm available in CardioViz3D.^{2} This method requires landmarks to be selected inside and outside the myocardium to fit an implicit elliptical surface. Second, the myocardium mask is meshed using CGAL software^{3}. This results in a tetrahedral mesh made up of around 65 500 linear tetrahedra. Both mechanical and electrophysiology models are highly related to the fibre directions. Realistic cardiac fibres are generated by synthetically varying the elevation angle across the myocardium wall. Regarding the left ventricle (including the septum as a part of it), the elevation angle is defined from +70° on the endocardium, to 0° at the middle of the endocardium and to −70° in the epicardium. Subsequently, the same method is applied to the right ventricle (figure 2).
2.1.2. Electrophysiology modelling of the heart: the Mitchell–Schaeffer model
As introduced previously, there exist several models to describe the ventricular AP. However, only phenomenological models should be considered for our objective of interactive simulation. Indeed, biophysical models characterize the electrophysiology at the cellular scale, which implies a high level of complexity and an important computational cost. Eikonal models are efficient algorithms but have difficulty simulating complex electrophysiology phenomena with nonphysiological parameters.
We chose the Mitchell–Schaeffer model because of its six parameters that can be physiologically identified. Moreover, it captures well the restitution properties of the APD compared with phenomenological models of similar complexity. The AP shape described by the Mitchell–Schaeffer model [9] is given in figure 3.
The Mitchell–Schaeffer model, derived from the Fenton–Karma model [18], has two variables: u, the transmembrane potential, and z, a secondary variable controlling the repolarization phase. The temporal evolution of these two variables is governed by the following equations: 2.1
The diffusion term is defined by a 3 × 3 anisotropic diffusion tensor , so that the planar conduction velocity (CV) in the fibre direction is 2.5 times greater than in the transverse plane (r = 1/(2.5)^{2}). d is the diffusion coefficient. The parameters τ_{in} and τ_{out} define the repolarization phase, whereas the constants τ_{open} and τ_{close} manage the gate opening or closing depending on the changeover voltage u_{gate}. The default values (describing the common AP) of these parameters are given by Mitchell & Schaeffer [9], and a simulation of an electrophysiological wave is shown in figure 4.
The term J_{stim}(t) is the stimulation current applied in the pacing area. Only a set of nodes corresponding to the apex of the endocardial surface of the left and right ventricles will include this stimulation current to initiate the depolarization wave (in sinus rhythm) in the ventricles. A secondary area can also be defined to model the stimulation induced by a catheter. For our simulations, we initiate the sinus stimulus during 0.1 ms so that 2.2
Based on the work of Coudiere and colleagues [19], extracellular potentials can nevertheless be estimated, because the authors have shown that it is possible to recover the extracellular potential from the transmembrane potential under a reasonable hypothesis. To do so, an additional elliptical equation has to be solved on the surface using an algebraic equation or an integral formulation. We can therefore easily estimate an extracellular potential using the monodomain Mitchell–Schaeffer model.
2.2. Numerical settings
We first focus on the spatial discretization method and the integration schemes under consideration. We subsequently detail a comparison study on the numerical settings of our simulation. Finally, we explain the way we handle the personalization of the Mitchell–Schaeffer model to make it patient specific.
2.2.1. Spatial discretization
As explained in §2.1.1, we compute our threedimensional volumetric mesh from MR images. The biventricular geometry is meshed with 65 500 linear tetrahedra, using CGAL. The Mitchell–Schaeffer model and the diffusion term are implemented, using the finite element method (FEM). The edge length for cardiac electrophysiology is usually less than 0.5 mm. However, we use larger tetrahedra (with an edge length of approx. 2 mm) in our simulations. Using larger elements means fewer tetrahedra within the mesh, thus ensuring higher computational performances. The influence of the element size will be studied in more detail in §2.2.3.
Recent work considered other spatial discretization methods. For instance, Rapaka et al. [20] presented their algorithm applying the lattice Boltzmann method to cardiac electrophysiology (LBMEP). In this work, the authors explained that traditional FEMs could not offer fast electrophysiology simulations. In our paper, we propose using the FEM model on a graphics processing unit (GPU) as a good alternative to compensate for the FEM computational cost.
2.2.2. Integration schemes
We now detail the integration schemes we considered. Ethier & Bourgault [21] presented the main schemes used for cardiac electrophysiology. We implemented the modified Crank–Nicholson/Adams–Bashforth (MCNAB), a secondorder, semiimplicit scheme. The diffusion term is implicitly integrated, whereas the ionic current term is explicitly computed. Electrophysiology equations can be generalized as follows: 2.3where g and h denote the ionic functions of the Mitchell–Schaeffer model, and f is the additional diffusion term (see equation (2.1)). The MCNAB scheme describes the integration as 2.4where M represents the operator obtained by integrating the term of mass density ρ; G and H denote the operator obtained by integration of the ionic term; and F denotes the operator obtained by integration of the diffusion term. In this last system, the index number n refers to the nth time step. Using the MCNAB, our simulation includes a conjugate gradient coupled with a Jacobi preconditioner in order to efficiently solve the linear system (Ax = b). Because the matrix A is diagonal dominant, the choice of this preconditioner is straightforward. The preconditioned matrix given by the Jacobi can be updated during the simulation in order to take into account conductivity changes (owing to thermoablation).
We also implemented a fully explicit scheme with a backward differentiation. In the following, we refer to it as ‘explicitBDF’. It can be solved without a linear solver. Using the notation from equation (2.3), it can be written as: 2.5
Maps of depolarization times obtained using both MCNAB and explicitBDF solvers are shown in figure 5a and figure 5b, respectively.
2.2.3. Numerical study
To make the best compromise between performance and accuracy, we study each of the numerical parameters of our simulation. The influence of the element size regarding the integration method is considered as well as the influence of the time step.
Element size. We first focus on a crucial feature regarding the performance context: the element size. Pathmanathan et al. [22] studied the relationship between the CV and the element size depending on the integration method. We reproduced the onedimensional planar propagation wave, based on the cell model of Mitchell–Schaeffer, on a simplified geometry: a threedimensional cuboid such as the one described by Pathmanathan et al. [22]. We implemented the lumped integration and the ionic current integration (ICI) that were described by Pathmanathan et al. [22]. The lumping method consists of summing all the coefficients of a line onto the diagonal. The ICI interpolates the nodal ionic current linearly on each element.
The results using the lumped integration and the ICI are presented in figure 6. We obtain the same trend as Pathmanathan et al. in their work. The ICI method is a more accurate technique than the lumped integration. However, the ICI integration method is more computationally demanding than the lumping method. We therefore decided to rely on a lumped integration of the ionic term of the Mitchell–Schaeffer model.
As presented in figure 6, the coarser the mesh, the slower the CV (with a constant diffusion coefficient). In our approach, we propose adapting the nodal diffusion coefficient d in order to fit the patientspecific map of depolarization times. The integration error owing to the lumping approximation would therefore be numerically compensated while benefiting from larger elements, i.e. better performances. This will be carried out during the electrophysiology personalization detailed in §2.2.4.
Time step. Depending on the integration scheme used, the time step is limited either owing to the diffusion term or owing to the ionic term, depending on whether the explicitBDF or the MCNAB scheme (implicit diffusion) is used. These stability conditions have already been tackled by Ethier & Bourgault [21] and by Coudière & Pierre [23]. It can be shown that:

— using the semiimplicit MCNAB scheme, the ionic current defines the limit of stability: where R denotes the ionic term;

— using the explicitBDF scheme, the diffusion term determines the limit of stability. The time step has to be dt ≤ min(dx/2 × CV), where dt is the characteristic length of the elements, and CV is the CV of the depolarization wave inside this tetrahedron. However, the limit owing to the ionic term is still active. This means that dt cannot exceed 0.286 ms whatever the scheme. Using the mesh described in §2.1.1, the stability limit using the explicitBDF scheme is dt < 0.13 ms.
These stability conditions are theoretical values. In simulations, it can be noted that the time step dt can, in some cases, be slightly higher than these theoretical limits.
Simulation benchmarks. To validate our simulations, several benchmarks have been performed using the explicitBDF solver and a lumped integration method. First of all, we want to show that using coarse elements (coarser than usual electrophysiology simulation) and increasing numerically the diffusion coefficients (to fit the depolarization times) leads to good performances while keeping the reference AP shape, CV and APD. Tests have been carried out on the same simplified geometry as described in section Element size (same as in reference [22]), using the model parameters presented by Mitchell & Schaeffer [9]. The beam is 2 cm long, and the section area is 4 mm^{2}. The simulation is initiated on a face at one extremity with a stimulus current (see equation (2.2)).
Three different simulations are performed:

— The first one is using fine finite elements and a fine time step: edge size dx = 0.1 mm, dt = 0.001 ms and with a reference diffusion coefficient d.

— The second one is using coarse finite elements and a very fine time step: dx = 1 mm, dt = 0.001 ms and using an adapted diffusion coefficient d_{adapt} to fit the same depolarization times as those obtained in the first configuration.

— The third one is using coarse finite elements and a large time step: dx = 1 mm, dt = 0.1 ms and using the adapted coefficient d_{adapt} as well.
The resulting AP shape is shown in figure 7. The curves show that the chosen time step (dt = 0.1 ms) seems to correctly capture the AP shape. Moreover, it appears from these results that the use of a coarser element (dx = 1 ms) does not affect the AP shape. Nevertheless, it is known that using coarser elements slows down the propagation phenomenon (figure 6). The diffusion coefficients need therefore to be adapted to compensate for this slowdown.
In table 1, we demonstrate that using a large time step (consistent with the stability limit) does not affect either the CV or the APD. The edge size used for these computations is dx = 1 mm, and the diffusion coefficient is adapted (d_{adapt}). From the results shown in table 1, it appears that the time step does not affect either the APD or the wave propagation (CV) using the semiimplicit MCNAB. However, using the explicitBDF scheme, the time step influences the CV. The personalization step will have to be performed with a determined time step if the explicitBDF scheme is used.
From now on, the following numerical settings are used:

— average edge size is around 2 mm,

— using the MCNAB scheme: time step dt = 0.2 ms, or using the explicitBDF scheme: time step dt = 0.1 ms.
2.2.4. Electrophysiology personalization
MR images have been acquired in a patient with a left bundle branch block (LBBB) and postmyocardial infarction scars. The occurrence of an LBBB implies a late activation (and therefore late contraction) of the left ventricle. Patterns of ventricular depolarization are therefore altered, which can result in prolongation of the depolarization, thus leading to ventricular desynchronization. Electrical asynchronicity often leads to an asynchronous contraction, which is detrimental to the cardiac ejection efficiency.
To simulate this patientspecific electrophysiology, we need to personalize the Mitchell–Schaeffer model from patientspecific data, such as contact and noncontact intracardiac electrical mapping. We use the method presented by Relan et al. [24] to obtain personalized conduction parameters (the diffusion coefficients d) as well as restitution parameters (the time constants τ_{in}, τ_{out}, τ_{open} and τ_{close}). By estimating the diffusion coefficients d, the discretization error can be compensated. This personalization allows us to obtain simulations with less than 10 ms of error with respect to the measured depolarization times, and to also recover the restitution parameters of the APD. The interested reader should refer to Relan et al. [24] for more detail.
2.3. Interactive electrophysiology using a graphics processing unit
As stressed by Rapaka et al. [20], the use of classic FEM is more computationally demanding than using any other nodal method such as the LBM electrophysiology. In this section, we detail our GPU implementation that leads to very interesting performances for cardiac electrophysiology simulations.
2.3.1. Basic graphics processing unit computing
A simulation on a GPU uses the graphics card of the computer to perform computations instead of using the central processing unit (CPU). However, efficient implementation on the GPU is very complex because every step of the computation has to be cleverly shared out among the memory. In our approach, we decided to rely on the CUDA toolkit (dedicated software for NVIDIA's GPU) to develop a GPU version of our electrophysiology model. The SOFA public framework is already interfaced with CUDA, and some codes for GPU computing are available. Similar computations could be run using OpenCL or other GPU models.
The GPU architecture consists of several multiprocessors able to carry out highly parallel tasks independently. To benefit from the GPU computational power, the complexity of implementation lies in defining an optimal distribution of independent threads while minimizing the memory access. This is far from being trivial and is totally different from a CPU parallel approach (such as domain decomposition or other strategies).
Single float precision is used because it is sufficiently accurate regarding the variables of the Mitchell–Schaeffer model. Moreover, Bartocci et al. [25] established the computational cost of the double precision: it has been estimated to be twice more computationally demanding than the single float precision.
2.3.2. Graphics processing unit optimizations
For the electrophysiology simulation, the main task is to implement efficiently the computations of the ionic current and the diffusion terms. Using the lumped integration method, the ionic term of the Mitchell–Schaeffer model is a nodal value (computed at each vertex separately). Therefore, this part of the computation is a highly parallel task. This implementation uses classic parallelization methods: each thread is dedicated to one vertex and computes the contribution of the ionic term for this vertex while ensuring a tiled access in memory.
Owing to strong neighbouring dependencies, the memory access during the computation of the diffusion term is very sensitive. The parallel implementation of the diffusion can lead to writing conflicts: two threads solving two adjacent edges may simultaneously write on the same point. This problem can be automatically handled using the atomic function available in the latest versions of CUDA (more than CUDA 2.0). However, these atomic functions are not the most efficient methods. An innovative and very efficient algorithm has been developed by Allard et al. [26] to tackle the writing conflict issue. We applied this algorithm, originally designed for FEM mechanical equations, to our electrophysiology model. This technique consists of dividing the diffusion computation into two separate kernels. Before the computation starts, the vertex–edge topology is saved. During the simulation, the contribution of the diffusion is first computed on all edges. Then, the contribution on each vertex is computed by accumulating the values computed on the edges using the topology information stored previously. For more efficiency, a parallel reduction can be implemented for this last accumulation step. Compared with the atomic functions proposed by NVIDIA, this method presents two main advantages: it is twice as fast and it does not depend on the CUDA version.
Solving the electrophysiology system results in several manipulations of large vectors. This step has a high potential for parallelism, because many small operations have to be conducted on a large number of nodes. We therefore focus on optimizations such as the summing vector or equality operations. However, these simple optimizations offer only a limited performance gain.
Minimizing the memory access latency ensures the most efficient simulation. Therefore, the communications between the CPU and GPU during the computation have to be removed. This assumes the implementation of the whole simulation using CUDA. The whole implementation allows a global acceleration of × 28 (using the explicitBDF scheme), which is a very satisfactory result. This will be detailed in §4.
3. Simulation of cardiac electromechanics
The myocardium is the heart muscle that contracts depending on the concentration of free calcium ions in the cells: this is called the myocardial contractility. In this section, we focus on simulating the mechanical behaviour of the heart and defining the electromechanical coupling. Existing models proposed in the past 20 years differ only in the choice of mechanical material, electrophysiology model or electromechanical coupling. Our approach is based on the BCS model [14] that was later improved by SainteMarie et al. [27] and Chapelle et al. [28].
This mechanical model does not include any stretchactivated ionic channels. This means that only control of the mechanics by the electrophysiology is possible and that there is no feedback from the mechanics to the electrophysiology. This coupling can be useful to study the mechanical response related to several stimulation patterns. Such an approach has been shown to be predictive in a limited number of heart failure cases by Sermesant et al. [4]. The mechanicsbased motion could also be used to generate clinical images from simulation.
3.1. The electromechanical model of Bestel–Clément–Sorine
The contraction of the myocardium depends on the cellular ion concentrations along preferred directions. The muscle is contracting along fibres which are bundles of myofibril (figure 8). These myofibrils can be detailed as a series of sarcomeres, which are complex structures, including:

— thin filaments (actin) and thick filaments (myosin) responsible for the contraction and relaxation,

— elastic filaments (titin) bounding sarcomeres and Zdiscs.
In the cardiac extracellular matrix, connecting tissues, mainly made up of collagen and elastin, are the surrounding fibres.
The study of SainteMarie et al. [27] demonstrated the good properties of the BCS model. This model, based on a multiscale physiological description, is consistent with the laws of thermodynamics.
The BCS model can be decomposed into two different parts:

— a passive isotropic viscohyperelastic component corresponding to the natural elasticity and friction of connecting tissues;

— an active component accounting for the contraction induced by the electrical impulses.
This has been implemented in SOFA [16] using a Euler implicit solver and is detailed in the following sections.
3.1.1. Passive nonlinear elasticity
The first component of the BCS approach is a passive hyperelastic material modelling the behaviour of the connecting tissues (extracellular matrix). In our simulation, we consider the Mooney–Rivlin model with isotropic properties. Our BCS formulation therefore assumes a transverse isotropy. The strain energy described by the Mooney–Rivlin material is given in equation (3.1), 3.1where C_{1}, C_{2} are material parameters that need to be determined, and K is the bulk modulus. and are the isochoric invariants of the right Cauchy deformation tensor, . The first invariant verifies , and the second invariant is , where , and J is the Jacobian .
Again, we implemented our passive hyperelastic material using the FEM. However, the strain energy is not computed using the classic Galerkin FEM formulation; instead, we used the multiplicative Jacobian energy decomposition (MJED) method presented by Marchesseau et al. [29].
This discretization of nonlinear hyperelastic materials on a linear tetrahedral mesh leads to faster stiffness matrix assembly than the classic Garlerkin formulation. It is based on the multiplicative decomposition of terms that depend on J from the terms that depends only on the invariants of the Cauchy deformation tensor. The MJED method proved to be approximately 2.7 times faster than the standard FEM on all hyperelastic materials, isotropic or anisotropic.
3.1.2. Active fibre contraction
As explained previously, the muscle fibres are bundles of myofibrils. In the BCS mechanical model, we describe the myocardial contraction with an active component that depends on the transmembrane AP, v.
As shown in figure 9, the active part can be divided into three elements. First, a viscosity element depending on parameter μ is in parallel with the contraction stress tensor, σ_{c}. The viscosity element accounts for the energy dissipation owing to friction inside the sarcomeres. The resulting stress can be written as .
Second, an elastic component with Young modulus E_{s} verifying is defined in series with the previous elements. This mimics the elastic behaviour of the filaments (titin) joining the sarcomeres to the Zdiscs. The addition of this component leads to a global passive stiffness corresponding to the two passive elements added in parallel, being therefore a transversally anisotropic material with a different stiffness along the fibre direction from that in its transverse direction. This approach corresponds to a classic Hill muscle model that has the advantage of correctly representing isotonic (for ejection and relaxation phases) and isometric (for isovolumetric phases) contraction.
The stress and strain are related by 3.2where is the projection of the Green–Lagrange tensor E on the fibre (). Based on the Huxley model [30] and on statistical mechanics, the mechanical behaviour of the nanoscopic myosin and actin filaments (figure 8) can be characterized by the differential equations (see equation (3.3)) at the macroscopic scale. These equations link the active stiffness k_{c} with the active stress τ_{c}, 3.3where α is a constant related to the crossbridge release owing to the high concentration rate, k_{0} is the maximum stiffness and σ_{0} is the maximum contraction. The reduction factor n_{0} enables the introduction of the fact that the maximum contraction depends on the fibre strain . This is called the Starling effect. Finally, the ‘potential’ v is a control variable derived from the electrophysiology model (Mitchell–Schaeffer model in our case), 3.4This variable v is related to a free calcium concentration and varies depending on the depolarization time T_{d} and the APD. k_{ATP} and k_{RS} are kinetic constants, respectively, describing the rate of myosin ATPase activity (responsible for contraction) and the rate of calcium pumped back into the sarcoplasmic reticulum (responsible for relaxation).
3.1.3. Haemodynamic model
The whole BCS model describes the complex mechanical behaviour of the heart tissue. However, the blood flow is not considered in the equation. A haemodynamic model needs to be defined in order to take into account the presence of blood inside the ventricles.
Chapelle et al. [28] proposed a constraint formulation to simulate the blood circulation. The blood circulation in each ventricle (and the associated constraint) can be written in four phases:

— Ventricular filling: the atrial pressure (P_{at}) is higher than the ventricular pressure (P_{v}) and drives the blood to flow from the atrium to the ventricle through the open mitral or triscupid valve: where K_{at} corresponds to a linear law.

— Isovolumetric contraction: the valves close and the ventricle contraction starts: where K_{iso} relaxes the usual isovolumetric constraint (q = 0).

— Ejection: semilunar (aortic and pulmonary) valves open under the ventricular pressure (P_{v}) (higher than the arterial pressure ) and the blood ejection occurs: where K_{ar} corresponds to a linear law.

— Isovolumetric relaxation: the valves close while the relaxation starts and the atrium slowly fills up with blood: .
This constraint is solved efficiently in SOFA using a prediction/correction algorithm, avoiding the addition of the unknown P_{v} as a state variable by projecting the corrected nodal velocities directly on the constrained space. The details of such an algorithm are explained by Marchesseau et al. [16].
3.2. Parameter calibration from medical images
In order to better fit the patientspecific heart motion, its deformation field is needed. To compute this motion field, we use the iLogDenons tool developed by Mansi et al. [31] using nonrigid registration. A dense nonlinear transformation is applied to our fourdimensional clinical data. Using elastic and incompressible regulizers in the registration, this tool even allows some components of the twist motion of the heart to be recovered. This motion field is then used to estimate the ventricular volumes in the calibration phase. Indeed, to globally estimate the mechanical parameters, a sensitivity study was first performed to detect which parameters are sensitive to which clinical observations.
The unscented transform is then applied and allows four parameters to be estimated in one iteration, minimizing the differences between observations extracted from the simulated volume and the measured volume [16]. This algorithm, represented in figure 10, builds a covariance matrix between the relevant parameters and the observations Z (in our case, the minimum of the left ventricle volume and the minimum and the maximum of its derivative) spread around some initial parameter set θ^{0}. The new parameters are then found to minimize the difference between the meansimulated observations and the measured observations Z^{obs} with 3.5
More parameters can be estimated when pressure curves are available, as shown by Marchesseau et al. [32].
4. Results
4.1. Realtime electrophysiology simulation
4.1.1. Performances
As shown in §2.3, the computation of electrophysiology can be accelerated with GPU computing. The hardware used in our GPU simulation is a GeForce GTX580 with 512 cores. The results obtained using GPU are compared with computation made on a CPU Intel Xeon Z3550. Table 2 summarizes the performances using either the explicitBDF solver or the MCNAB solver. Performance is evaluated with the ‘realtime ratio’ (r_{RT}), which can be expressed as the ratio of the duration of the computation against the time that is simulated.
It can be noted from table 2 that the performance gain is less important using the MCNAB than using the explicitBDF scheme. Solvers (including preconditioner) used with the MCNAB require a lot of matrix–vector multiplications that have a high computing cost using the GPU. We can see from table 2 that using the explicitBDF solver allows realtime performances to be reached. This means one cardiac cycle () can be computed within less than 0.97 s.
4.1.2. User interaction
In the context of a training system for junior electrophysiologists, it is important that the user can interact with the electrophysiology simulation. In our framework, stimulations through pacing leads can be applied anywhere on the endocardium. The user can also define the power and the duration of the induced current J_{stim}. An example of stimulation at the apex of the heart is presented in figure 11.
The electrophysiology simulation can also handle thermoablation. As soon as cardiac cells are heated by the radiofrequency probe, their conductivity (diffusion coefficient) is then decreased to a null value. This means that, after ablation, we considered that the targeted cells no longer conduct an electrical current. However, our simulation would allow an intermediate ablation level by defining an intermediate conductivity. It is also possible to simulate the catheterized measurement of extracellular potentials from APs using the work of Coudière and colleagues [19].
4.2. Coupled electromechanical model
4.2.2. Interactive mechanics
Simulation of the whole mechanical model (BCS model with haemodynamical constraint) is sequentially computed. In our approach, cardiac mechanical contraction is solely driven by electrophysiology, thus ignoring a mechanoelectric feedback.
The dynamic equations of motion are solved here, on the same tetrahedral mesh as used for the electrophysiological model. They are discretized using a Euler implicit scheme and then solved using the conjugated gradient algorithm. First, we solve the electrophysiology on a static mesh that gives the APs for each node. Thresholding these APs gives depolarization times that are used as input for the mechanical simulation. The computation time of the cardiac mechanics during one cardiac cycle requires (for our 65 500 linear tetrahedral mesh) about 10 min. Nevertheless, the simulation being tractable, which allows us to consider efficient personalization strategies. Screenshots of the simulation are shown in figure 12.
4.2.2. Realistic simulation
By using personalized realtime electrophysiology with a calibrated mechanical model, realistic simulation of the cardiac cycle can be achieved. The late activation (i.e. late contraction) owing to LBBB can be observed, for instance, in figure 12. Our approach also allows us to produce realistic pressure–volume loops, as plotted in figure 13.
5. Discussion and conclusion
In this paper, we first presented an interactive and realtime simulation of cardiac electrophysiology. Based on electrophysiology simulations (depolarization times and APD), a tractable mechanical model of the myocardium has been developed. The results given in §4 confirm the ability of our framework to efficiently simulate patientspecific electrophysiology and mechanics. For instance, in the case of LBBB, late activation (or contraction) of the left ventricle is simulated leading to an asynchronous myocardial motion.
Future work will take advantage of the fast electrophysiology simulation in order to efficiently personalize the Mitchell–Schaeffer electrophysiology directly in the SOFA framework. Further optimization can be achieved by improving the numerical integration of the equations as well as by using new GPU hardware. Moreover, nonrealtime electrophysiology simulations could be considered using finer elements (dx ∼ 0.1 mm) in order to reproduce some more complex behaviours of patient electrophysiology. Such simulations would be interesting because they should remain fast (with good performance), thanks to our efficient GPU implementation. Nevertheless, nonrealtime simulations are incompatible with a training simulator.
Regarding the mechanical model of the heart, other haemodynamic constraint formulations could be envisioned, such as the CircAdapt constraint of Lumens et al. [33] that links the blood flow and pressure between the right and left heart. Also, a partial GPU implementation of the BCS electromechanical model could be achieved to significantly decrease the computation times of the mechanical part. Finally, the development of a training system for CRT and thermoablation coupling realtime endovascular navigation and electrophysiology simulation is underway to overcome the limitations of current platforms.
Acknowledgements
This work is partially funded by the European project euHeart FP720072013224495.
Footnotes
One contribution of 25 to a Theme Issue ‘The virtual physiological human: integrative approaches to computational biomedicine’.
↵1 More information about SOFA is available at http://www.sofaframework.org.
↵2 CardioViz3D is an open source software including processing, simulation and visualization tools for cardiac images. It can be downloaded from http://wwwsop.inria.fr/asclepios/software/CardioViz3D/.
↵3 The Computations Geometry Algorithm Library (CGAL) is available at http://www.cgal.org.
 Received November 12, 2012.
 Accepted January 8, 2013.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.