Royal Society Publishing

Coupled personalization of cardiac electrophysiology models for prediction of ischaemic ventricular tachycardia

Jatin Relan, Phani Chinchapatnam, Maxime Sermesant, Kawal Rhode, Matt Ginks, Hervé Delingette, C. Aldo Rinaldi, Reza Razavi, Nicholas Ayache

Abstract

In order to translate the important progress in cardiac electrophysiology modelling of the last decades into clinical applications, there is a requirement to make macroscopic models that can be used for the planning and performance of the clinical procedures. This requires model personalization, i.e. estimation of patient-specific model parameters and computations compatible with clinical constraints. Simplified macroscopic models can allow a rapid estimation of the tissue conductivity, but are often unreliable to predict arrhythmias. Conversely, complex biophysical models are more complete and have mechanisms of arrhythmogenesis and arrhythmia sustainibility, but are computationally expensive and their predictions at the organ scale still have to be validated. We present a coupled personalization framework that combines the power of the two kinds of models while keeping the computational complexity tractable. A simple eikonal model is used to estimate the conductivity parameters, which are then used to set the parameters of a biophysical model, the Mitchell–Schaeffer (MS) model. Additional parameters related to action potential duration restitution curves for the tissue are further estimated for the MS model. This framework is applied to a clinical dataset derived from a hybrid X-ray/magnetic resonance imaging and non-contact mapping procedure on a patient with heart failure. This personalized MS model is then used to perform an in silico simulation of a ventricular tachycardia (VT) stimulation protocol to predict the induction of VT. This proof of concept opens up possibilities of using VT induction modelling in order to both assess the risk of VT for a given patient and also to plan a potential subsequent radio-frequency ablation strategy to treat VT.

1. Introduction

Cardiac arrhythmias including ventricular tachycardia (VT) are increasingly being treated by radio-frequency (RF) ablation. These procedures can be very effective but still have unsatisfactory success rates widely ranging from 50–90%, with a 20–40% late recurrence rate, owing to a lack of clinical consensus on the optimum RF ablation strategy [1]. There is a need for substantial guidance in locating the optimum ablation strategy [2].

This guidance could be provided by personalized in silico cardiac electrophysiology models, as such models may allow different ablation strategies to be tested. A personalized model incorporates estimation of patient-specific parameters, which best fit the clinical data. Such a step is necessary to reveal hidden properties of the tissue that are used to predict the behaviour under different pacing conditions.

There is a large variety of cardiac electrophysiology models for myocyte action potential developed at cellular and sub-cellular scales [37]. Cardiac tissue and whole-heart electrophysiological computations of these models are based on the principles of reaction–diffusion systems [8,9]. According to the reaction term computation, these models can be broadly categorized as biophysical models (BMs), phenomenological models (PMs) and generic models (GMs). BMs [5,6] model ionic currents and are the most complete and complex but are less suitable for parameter estimation from clinical data owing to a high computational cost and to the lack of observability of their parameters. PMs [10,11] are based on PDEs and are of intermediate complexity level and less computationally expensive. GM [12,13] represent simplified action potentials and are the least complex. Simple eikonal (EK) models [1416] model the action potential propagation in the cardiac tissue without modelling the action potential itself. They can be very fast to compute [17], but less reliable in arrhythmia predictions due to the complexity of both the refractoriness and the curvature of the wavefront.

Computational modelling of cardiac arrhythmogenesis and arrhythmia maintenance using such models has made a significant contribution to the understanding of the underlying mechanisms [1823]. These studies have shown a host of factors involved in the onset of arrhythmia with wave fragmentation and spiral wave breakups, which include realistic ventricular geometry [24], heterogeneity in repolarization [25] and action potential duration (APD) restitution [26,27] and conduction velocity (CV) restitution [28]. A combined clinical study and synthetic modelling of APD restitution was shown in Nash et al. [29]. In this paper, we study these properties for a clinical dataset and evaluate its role in ischaemic VT induction.

To introduce models directly into clinical practice, the ideal requirements are low computational complexity, fast estimation of parameters (quick personalization) and reliable predictions. These attributes cannot be found in one single model, thus here we present a novel approach, wherein we combine two models to obtain these attributes and apply them to a clinical dataset. We use a coupled personalization framework, which is fast and combines the benefits of an EK model [30] with those of a simplified BM, the Mitchell–Schaeffer (MS) model [31]. The fast three-dimensional EK model is used to estimate the tissue conductivity parameter over the ventricles derived from non-contact mapping of the endocardial surface potential using an adaptive iterative algorithm. This is then used to set the conductivity parameter of the three-dimensional MS model. Additional parameters related to APD restitution properties of the tissue are then estimated locally using directly the three-dimensional MS model and the measured endocardial surface potential.

This framework is applied to a clinical dataset from a patient with heart failure and myocardial scar on magnetic resonance imaging (MRI) scanning using electrophysiological data from a non-contact mapping study performed in a hybrid X-ray/magnetic resonance (XMR) suite [32]. The ventricles were mapped with a statistical atlas for cardiac fibres [33] (figure 1). The resulting personalized three-dimensional MS model is then used to simulate a clinical VT-stimulation (VT-stim) procedure to show a potential application of VT induction modelling. Figure 2 shows the framework of the coupled personalization method and VT induction modelling used in this paper.

Figure 1.

(a) MR derived segmented mesh with scars (in red). (b) Fibre orientations based on a statistical atlas.

Figure 2.

Flowchart describing the outline of this paper.

2. Clinical context

In order to evaluate the inducibility of VT in patients, the clinical procedure involves stimulation with an EP catheter usually in the right ventricle apex (RVA) at different cycle lengths. This type of protocol is used to test if re-entrant VT can be induced by such pacing in patients at risk of VT. Such studies may be useful in predicting the risk of VT for an individual patient but provide limited information on which to base a potential ablation strategy of re-entrant VT circuits.

Our aim is to create a personalized electrophysiological model of a given patient to which a virtual VT-stimulation (VT-stim) procedure can be applied. Moreover, a virtual RF ablation procedure can then be applied to the model in order to test potential ablation strategies.

In order to validate this approach, we have used an extensive clinical dataset, derived from mapping of the left ventricle (LV) endocardium. Such mapping is not routinely performed for a VT-stim procedure; however, it is sometimes used to guide RF ablation. This clinical dataset was obtained from an electrophysiology study performed in a hybrid X-ray and magnetic resonance (XMR) environment. The electrical measurements obtained using an ensite system (Saint Jude Medical) were registered to the patient anatomy using XMR registration [32] (figure 3b).

Figure 3.

(a) T-wave polarity map on ensite LV surface, with positive (red), negative (green), bi-phasic (blue) and undetected (black) T-waves. (b) Repolarization times projection from ensite LV surface to MR LV endocardium after XMR registration. Black, undetected repolarization times.

2.1. Depolarization and repolarization times extraction

The electrical data were collected with high pass filter settings for prominent QRS detection and with low pass filter for T-wave detection.

The depolarization times (DTs) were detected within the QRS window set from surface electrocardiograph (ECG; figure 4a), and were derived from the zero crossings of the Laplacian of the measured unipolar electrograms V (figures 4b and 5a). The surface Laplacian of electrograms gave reasonable results even in the case of multiple deflections and also allows detection of local activation without interference from far-field activity [34].

Figure 4.

(a) Measured surface ECG I, with QRS window (for DT extraction) and ST window (for repolarization time extraction). Ensite surface unipolar electrograms (b) with high-frequency band-pass filter for detection (black dots) of DTs, and (c) with low-frequency band-pass filter for detection of repolarization times, from positive (red), negative (green) and biphasic (blue) T-waves. Few electrograms had indiscernible T-waves (black).

Figure 5.

(a) Comparison of the measured DT isochrones on the LV surface only with model simulated DT isochrones on the whole heart after personalization (b) shows the same for measured (LV surface only) and model simulated (whole heart) APD maps. Measurements are for baseline.

The repolarization times were detected within the ST window in surface ECG (figure 4a) and were derived using the ‘alternative method’ when compared with the standard Wyatt method, because (i) most of the positive T-wave electrograms had an indiscernible steep upstroke for repolarization time detection with the Wyatt method, (ii) a closer correlation was obtained between activation recovery interval (ARI) and monophasic action potential (MAP) duration with the alternative method, as discussed in Yue et al. [26], and (iii) the difference in APD extraction from the two methods had only a minimal influence on the restitution slopes and spatial distributions [29]. The measured T-wave polarity maps had positive T-waves for early repolarizing sites and negative for late, as in agreement with Potse et al. [35] (figure 3).

The alternative method has repolarization times derived from dV/dtmax for the negative T-wave, at the dV/dtmin for the positive T-wave, and the mean time between dV/dtmax and dV/dtmin for the biphasic T-waves (figures 3, 4c, 5c and 6).

The data were collected during intrinsic sinus rhythm and atrial pacing mode at 100 beats per minute (b.p.m.; figure 6). Myocardial scar was segmented manually from the delayed enhancement MR image.

Figure 6.

Steady-state APD values measured from the ensite data and projected on the LV surface, for (a) baseline and (b) pacing in the atria (at 100 b.p.m.). Black represents scars.

The patient was a 60-year-old woman with heart failure and New York Heart Association (NYHA) class III symptoms. The patient had a dilated cardiomyopathy with subcritical disease on coronary angiography, although cardiac MRI showed subendocardial postero-lateral scar in the LV. The left ventricular ejection fraction was 25 per cent on maximal tolerated heart failure medication. The surface ECG demonstrated significant conduction disease with left bundle branch block and a QRS duration of 154 ms (normal QRS is less than 120 ms). Echocardiography, including tissue doppler, confirmed significant mechanical dysynchrony in keeping with the ECG findings.

This patient was selected for cardiac resynchronization therapy (CRT). When implanting a CRT device, there is a choice between a standard device and one also integrating a defibrillator. Therefore, evaluating the risks of VT in such CRT patients is important.

3. Cardiac electrophysiology models

3.1. Eikonal model

The EK model simulates the propagation of the depolarization wave in cardiac tissue, ignoring the repolarization phase. The EK model is governed by the EK diffusion equation [15,36] and is solved using the fast marching method (FMM). It can be written asEmbedded Image 3.1where the superscript t denotes transpose, c0 is a dimensionless constant (=2.5), and τ(x) is the cell membrane time constant (=0.003 s). dEK is the square of the tissue space constant along the fibre and is related to the specific conductivity of the tissue in the fibre direction and has units of m2. The anisotropy is incorporated in the diffusion tensor M = diag (1,ρ, ρ), with ρ the anisotropy ratio between longitudinal and transverse diffusion. We use ρ = 1/2.52 in order to have CV 2.5 times faster in the fibre direction [14]. The nonlinear equation (3.1) is solved using a fixed point iterative method combined with a very fast EK solver as explained in detail in [17,30].

3.2. Simplified biophysical model (Mitchell–Schaeffer model)

The MS model [31] is a two-variable simplified BM derived from the three-variable Fenton Karma (FK) ionic model [37]. It models the transmembrane potential as the sum of a passive diffusive current and several active reactive currents including combination of all inward (primarily Na+ and Ca2+) and outward (primarily K+) phenomenological ionic currents. The MS model is described by the following system of partial differential equations,Embedded Image 3.2where u is a normalized transmembrane potential variable, and z is a gating variable, which makes the currents gate open and close, thus depicting the depolarization and repolarization phase. Jin = (zu2 (1 − u))/τin represents the the inward currents which raise the action potential voltage and Jout = −u/τout represents the outward currents that decrease the action potential voltage describing repolarization. Jstim is the stimulation current, at the pacing location. The literature values for reaction term parameters are given in Mitchell & Schaeffer [31]. τin, τout, τopen and τclose has units of s. The diffusion term in the model is also controlled by the diffusion tensor M. In the longitudinal direction of the fibre, this pseudo-conductivity is set to dMS, which is one of the parameters we adjust, and to dMS · ρ in the transverse directions. dMS has units of s−1. The electrophysiology model is solved spatially using the P1 finite element method (FEM), and in time using an semi-implicit scheme as modified Crank–Nicolson/Adams–Bashforth (MCNAB) scheme, which is evaluated in terms of accuracy, stability and computational time [38]. The parameter values and simulation details are given in tables 1 and 2.

View this table:
Table 1.

Electrophysiology model parameters.

View this table:
Table 2.

Model simulation and personalization specifications.

Scars were modelled with zero conductivity in the ischaemic zones. While the grey zones and the regions between scars (isthmus) had conductivity, APD restitution curves (RCs) estimated from the data, as shown in §5a.

4. Coupled personalization method

4.1. Apparent conductivity parameter estimation

Cardiac tissue conductivity is a crucial feature for the detection of conduction pathologies. The apparent conductivity (AC) of the tissue can be measured by the parameter dEK in the EK model.

4.1.1. Left ventricle endocardial values

AC is initially estimated on the endocardial surface as a global value using a simple bisection method which matches the average conduction velocity of the measured DT isochrones to the simulated ones. Using it as an initial guess, an adaptive multi-level domain decomposition algorithm is used, which minimizes the mean-squared difference of the simulated and measured DT isochrones at each level using a Brent's Optimization Algorithm presented in Chinchapatnam et al. [17]. Owing to the absence of transmural electrical propagation information, we assume no variation across the LV myocardium (excluding LV endocardium and scars) and hence, we prescribe a single value for the myocardial tissue across the LV wall.

4.1.2. Left and right ventricle myocardial values

The AC values for RV endocardium and RV myocardial mass are set at 5 and 0.64 mm2 (from the literature [14]). The LV myocardial AC value is estimated by one-dimensional minimization of the following cost function (mean-squared difference of simulated and measured isochrones at endocardium + squared difference of simulated and measured QRS duration). The simulated QRS duration is calculated as the difference between the maximum and the minimum DTs in the biventricular mesh and the measured QRS duration is estimated from the surface ECG.

4.2. Coupling of eikonal and Mitchell–Schaeffer model parameters

The AC parameter for EK model dEK (equation (3.1)) is a scale for the diffusion speed of the depolarization wavefront in the tissue. The model conduction velocity (CV) is related to dEK (figure 7) asEmbedded Image 4.1where the constants αEK and βEK are introduced to take into account the discretization errors (in particular of the curvature) in three dimensions.

Figure 7.

Parameter dEK and dMS relationship with simulated conduction velocity CV for both models. Blue solid line, EK model CV curve—fitted; blue open squares, EK model CV data—simulated; black solid line, MS model CV curve—fitted; black open circles, MS model CV data—simulated.

The corresponding conductivity parameter for MS model, dMS is also a scale for the wave diffusion speed in the tissue. The model CV here is related to dMS (figure 7) as,Embedded Image 4.2where the constants αMS and βMS are introduced for the same reasons as for the EK model, while τin is kept as a constant. The estimated AC parameter dEK can then be used to estimate the parameter dMS. The parameter dEK gives EK model CV, cEK, which is similar to the actual measured data CV (cmsd) after the parameter estimation step. Thus to have MS model CV (cMS) similar to the measured data, it has to be similar to EK model CV (cEK). The constants αEK and βEK represent numerical curvature, diffusion and discretization errors for the EK model based on FMM. They are different from the constants αMS and βMS, which are diffusion and discretization errors based on FEM. These constants are determined in three dimensions for the ventricular mesh. We performed several simulations with various dEK and dMS values and noted the corresponding cEK and cMS values. Then, we fitted the analytical curves given in equations (4.1) and (4.2) in a least square sense and determined the constants. The constants estimated are αEK = 802.25, βEK = −268.54, αMS = 995.87 and βMS = −554.38. Thus from equations (4.1) and (4.2), we have Embedded Image and Embedded Image. For a given CV, dMS is slightly different from dEK. This may be due to the fact that d has different units, as for the EK model we model DT and for the MS model we model transmembrane potential. However, they are both scales for diffusion, thus such a coupling is performed. Then, the personalized dMS values are computed from the corresponding estimated dEK values using equation (4.3) based on the condition that cmsd = CV = cEK = cMS after personalization. Thus from equations (4.1) and (4.2), we haveEmbedded Image 4.3

4.3. Parameter estimation for action potential duration restitution

APD restitution is an electrophysiological property of the cardiac tissue and defines the adaptation of APD as a function of the heart rate. Its slope has a heterogeneous spatial distribution, which can play a crucial role in arrhythmogenesis [9,22,29]. The APD RC defines the relationship between the next cycle APD and the diastolic interval (DI) of the previous cycle. The slope of these RCs is controlled by τopen and depicts the APD heterogeneity present at multiple heart rates. APD RC for MS model is explicitly derived as [31],Embedded Image 4.4where hmin = 4(τin/τout) and n is the cycle number. The maximum value of APD is also explicitly derived as,Embedded Image 4.5

4.3.1. Left ventricle endocardial values

We had recordings for the paced mode with 100 b.p.m. and a sinus rhythm rate (baseline). Therefore, for the estimation of the RC slope from APD-DI values at two rates, we assume that APDmax (asymptotic value of APD RC) should be approximately equal to the normal sinus rhythm APD, and the slope value is adjusted with paced mode APD. From equation (4.4) and (4.5), we can observe that τclose and hmin control both APD RC and APDmax, hence we estimate the parameters minimizing the error on them jointly. The cost function minimized isEmbedded Image 4.6with N the total number of pacing rates, i the vertex having data (LV surface only), θ = [τclose, τopen], APDSRmsd as measured by sinus rhythm APD and DImsd measured from the data as DImsd = 1/f − APDmsd, where f is the heart rate (in hertz).

Only APD restitution is estimated as we have two pacing frequencies, but additional data would allow to also estimate CV restitution, as demonstrated by experimental data in Relan et al. [39,40]. Parameter hmin is not estimated here but kept to the literature value [31] as it also controls the CV, unlike τclose and τopen, thus disturbing the CV adjustments done before. The parameter optimization method used here is a nonlinear constrained Active-Set Algorithm, with constraints on τclose and τopen to be in the range of literature values [31]. Figure 8c shows the fit of RC to data.

Figure 8.

Estimated parameters: (a) conduction velocity estimated from AC maps, (b) APD parameter τclose, lower τclose values correspond to lower measured APD (white ellipse), (c,d) heterogeneous APD RCs and APD RC parameter τopen maps, respectively, low τopen values (red) correspond to steep RC slopes and high values (blue) correspond to flat RC slopes. Heterogeneous minimum DI values for the RCs are also shown in (c).

Minimum DI can also be computed explicitly from the estimated parameter values (figure 8c) as described in Mitchell & Schaeffer [31],Embedded Image 4.7whereEmbedded Image with Jstim as defined in equation (3.2).

4.3.2. Left and right ventricle myocardial values

For RV, we fix one value measured from the QT interval given through the surface ECG. To have a smooth gradation of APD restitution from epicardium to endocardium, we diffuse the τclose and τopen values spatially in the LV myocardium from endocardium to epicardium as in Keldermann et al. [41] (figure 8d).

Steady-state APDs for a pacing frequency f could then be estimated from the intersection of the line DIn = 1/f − APDn+1 with the personalized RCs (figure 9). The model personalization times are given in table 2.

Figure 9.

Computation of local steady-state APDs for different cycle lengths from the estimated RCs. Red to blue colours represent steep to flat slopes.

5. Results

5.1. Parameter estimation

The AC parameters estimated using the EK model (figure 8a) show high conduction areas on parts of the endocardium, potentially depicting the Purkinje network extremities, and a conduction block near the scar. The coupled MS model conductivity parameters were then estimated from this AC map. The mean absolute error on simulated DTs after personalization was 7.1 ms for the EK model and 18.5 ms for the MS model (≈10–14% of depolarization duration 131 ms; figures 5 and 10a). The mean absolute error on APD was 8.71 ms (≈2% of APD 300 ms), showing a good fit as well (figures 5 and 10b).

Figure 10.

Maps of absolute (a) DT error and (b) APD error between simulated and measured isochrones after personalization.

Figure 8b shows the heterogeneity of the steady-state APD in terms of the estimated parameter τclose, as it is shorter on the free wall of the LV compared with the septum (white ellipse). Also there is a longer APD compared with the neighbours near the scar (grey zones) and the region between the two scars (isthmus; figure 8b (black ellipse)). For APD restitution, the mean absolute error after fitting the curves was 1.13 ms, showing a good fit (figure 8c). The region around the scars and the isthmus were more heterogeneous for the RC slope parameter τopen than the rest of the LV (figure 8d, black ellipses). In figure 8d, red (low values of τopen) corresponds to steep slopes and blue (high values of τopen) corresponds to flat slopes, as shown in figure 8c. The colour bar is also applicable to the RCs. A smooth apex-to-base gradient for APD RC can be observed in figure 8d. Isthmus and grey zones were seen to have high APD values with more heterogeneous RC slopes and conduction. All of these factors along with the scar geometry played a crucial role in induction of ischaemic VT as explained later.

5.2. Assessment of heterogeneity maps

Such maps reveal areas with a risk of VT (black ellipse in figure 11). But this is also highly dependent on the pacing location, thus in order to really assess this risk we simulate VT-stim procedures with various pacing locations.

Figure 11.

Steady-state values of APD for different pacing rates. Black ellipse highlights the changes in heterogeneity near the scars and isthmus.

In order to better interpret these parameter maps, we estimated the steady-state APD maps for different pacing frequencies (figure 11). We can observe that this map is quite homogeneous for slow pacing frequency (60 b.p.m.), but its heterogeneity then increases with pacing rate (100, 150 b.p.m.). Eventually it reaches a kind of plateau where it is then quite homogeneous (200 b.p.m.).

5.3. Ischaemic ventricular tachycardia stimulation

Programmed ventricular stimulation is a clinical protocol and consists of a number of extra stimuli introduced at two ventricular sites (RVA and RV-outflow tract (RVOT)), using various cycle lengths (CL), with varying coupling intervals. This protocol is tested directly on the patient, without any planning, to collect information about the VT and to plan the RF ablation lines. We illustrate in figure 12 how such pacing can develop VT. It may be time consuming or fail when VT is not inducible and recurrent. Here we followed a conventional VT-stim protocol [42] with an RVA pacing site, two extra stimuli and a shortest coupling interval of 200 ms at 600 ms pacing cycle length. The induced VT specifications are given in table 3.

View this table:
Table 3.

Induced ventricular tachycardia specifications.

Figure 12.

DT isochrones for simulated S1-S2 VT-stim protocol. Personalized anatomy and synthetic electrophysiology parameter set were used here with steep RC slopes and low CV in isthmus compared to flat RC slopes and high CV in the rest of myocardium. S1 stimulus shows a normal propagation and S2 shows a unidirectional block created in the isthmus due to APD heterogeneity. Bottom row shows DT isochrones for induced monomorphic VT. Figure adapted from Relan et al. [43]

5.3.1. Synthetic simulation

The first simulation of this protocol in silico was performed on the patient's ventricular and scar geometry, using a synthetic model parameter set as shown in figure 12, and detailed in Relan et al. [43].

5.3.2. Personalized simulation

We then used the personalized three-dimensional MS model of our patient data to simulate this protocol. The VT-stim simulation using RVA and RVOT on this patient did not induce any VT. This is in agreement with the clinical data on this patient who did not have any VT episodes. However, we could induce VT using various other pacing sites in LV, as explained next.

5.4. Towards ventricular tachycardia risk maps

Personalized models offer much more flexibility in the VT-stim procedure, as the model can be paced from any location and at any pacing rate, which would not be possible in practice. We used this by running a series of VT-stim procedures from many different points of the heart on a cluster of computers1. This provides insights on which points of the heart could be a source of VT if paced [44]. We then obtained sustained VTs when pacing from particular locations (figure 13),2 due to the surrounding heterogeneity in restitution and conductivity (figure 8).

Figure 13.

DT isochrones for simulated VT-stim protocol for over-drive pacing using the personalized electrophysiology parameters, with pacing near the scars (shown by arrow) at 400 ms cycle length ( 150 b.p.m.). (a) shows normal propagation (b, c) show increase in heterogeneity in conduction near the isthmus after seven cycles due to APD and refractoriness heterogeneity (black ellipse) and (d) no pacing and induced VT of 250 ms cycle length (≈ 240 b.p.m.) with an exit point (shown by a star).

Figure 14.

Induced VT in the personalized model. At the seventh cycle of the pacing overdrive, some areas in the isthmus have the previous S1 wave-back touching the wave-front of the next S1 wave, thus creating an uni-directional block. A monomorphic sustained VT then develops.

We now need to obtain clinical data on activation maps of sustained VTs in order to validate our VT simulations [44,42]. Such maps are difficult to acquire as routine VT-stim procedures do not include electro-anatomical mapping. We are currently planning a study where such data will be acquired. From these VT-stim simulations we can build a VT risk map for a given patient. As an ectopic focus can be generated by ischaemic cells [45,46], thus from any location in the myocardium, building such risk maps could have an important clinical impact when taking decisions regarding defibrillator implantation and RF ablation.

6. Discussion

6.1. Data limitations

As only LV endocardial mapping data were used in this paper, the model personalization had several limitations, a few of which are (i) lack of estimation of local RV and transmural spatial distributions of conductivity and restitution properties. A global estimation of these properties was done using body surface ECG waveforms as explained in §4a,c; (ii) only two heart rates were used for fitting APD RCs. Thus only two model parameters were estimated, and an assumption was made by considering the sinus rhythm APD as a constraint on the asymptotic value of APD RCs controlled by the parameter τclose, while the paced mode APD was used to adjust the RC slope controlled by the parameter τopen. However, more measurements for frequencies in the slope region (region between minimum DI and maximum value of APD on APD RC) could have depicted the RC slope more accurately; (iii) lack of minimum DI measurement with pacing periods reaching up to the refractory period of the myocardium. As these values can play a crucial role in the initiation and sustainability of arrhythmias. However, the model had its own simulated minimum DI derived from the estimated parameters as shown in §4c; and (iv) usage of non-contact mapping (NCM) data and EP and MR fusion errors. Although NCM data do have an advantage of measuring temporal EP data with more spatial acquisition (surface) than the contact mapping data (point), the NCM data can be challenging for local depolarization and repolarization time estimations. Also uncertainty on the data can be added due to the difficult registration between the ensite LV surface and the MR-derived LV surface (figure 3).

6.2. Model simplifications

In order to have a clinically relevant model for personalization and VT risk assessment, the MS model used had several simplifications, some of which are (i) no actual Purkinje network modelling: exact locations of these Purkinje network extremities are ambiguous and inextractable from the patient's data, although they could play a crucial role in arrhythmogenesis [46,47]. However, a high conductivity endocardial region was obtained, which may be inferred as depicting the underlying Purkinje network with personalization (§4a and figure 8a). A local estimation of endocardial restitution properties (§4c and figure 8d) also helped potentially depict the abnormalities in the Purkinje network around scars, leading to arrhythmia generation. (ii) Use of an atlas-based cardiac fibre model: extraction of true in vivo cardiac fibre orientations is the subject of ongoing research and including them would give more accuracy to the VT-stim predictions and inducibility maps. Finally, orthotropic anisotropy could change the model behaviour, but acquiring patient-specific data on cardiac laminar sheets seems even more challenging.

7. Conclusion

The proposed approach of coupled model personalization for fast estimation of hidden parameters such as conductivity and APD restitution could enable the clinical use of cardiac electrophysiology models in the future. The parameter estimation algorithm is used on clinical interventional data and the obtained results are very encouraging. The estimated conductivity and APD restitution parameters are able to distinguish between the healthy areas and the pathological ones (scar and isthmus). A clinical VT-stim protocol was simulated on the personalized MS model in order to assess the risk of VT and fibrillation for a patient-specific case. This opens up possibilities of evaluating the role of patient-specific models in the clinical setting to provide aid in treatment and planning of RF ablation. In the future, we need to validate the personalized model predictions with mapping data for arrhythmias and integrate the simulation of RF ablation into a first cohort of clinical cases.

Acknowledgements

This research received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 224495 (euHeart project).

Appendix

The MS model was personalized on a volumetric tetrahedral mesh of a four valve biventricular anatomy with a mean edge length (MEL) of 1.0 mm and 65 547 number of tetrahedrons. The model parameters, simulation, personalization and VT induction specifications are detailed in tables 1,2 and 3, respectively.

Footnotes

  • One contribution of 17 to a Theme Issue ‘The virtual physiological human’.

  • 1 64 nodes of IBM e325 dual-opterons 246, 2 GHz.

  • 2 A video on personalized induced VT simulation is available as electronic supplementary material. A few snapshots of the same are presented in figure 14 demonstrating the re-entry wave.

  • Received November 22, 2010.
  • Accepted March 3, 2011.

References

View Abstract