Fast personalized electrophysiological models from computed tomography images for ventricular tachycardia ablation planning

Abstract

Aims

Clinical application of patient-specific cardiac computer models requires fast and robust processing pipelines that can be seamlessly integrated into clinical workflows. We aim at building such a pipeline from computed tomography (CT) images to personalized cardiac electrophysiology (EP) model. The simulation output could be useful in the context of post-infarct ventricular tachycardia (VT) radiofrequency ablation (RFA) planning for pre-operative targets prediction.

Methods and results

The support for model personalization is a patient-specific virtual three-dimensional heart obtained from CT images. Here, the scar is identified as thinning of the myocardial wall on automatically computed thickness maps. We then use an Eikonal model of wave front propagation with reduced velocity in the damaged areas. An image-based vessel enhancement algorithm can automatically identify VT isthmuses. The personalized model is used for virtual pacing. We obtained a very fast pipeline that enables simulations in only a few minutes. It is fully automated starting from the semi-automated image segmentation phase. The computational time frame is compatible with the construction of a virtual pacing tool. In this tool, onset points and an optional directional block could be interactively selected. The directional block is a simple way to model tissue refractoriness. Output activation maps are compared with EP data acquired pre-operatively. We show that this framework allows the reproduction of recorded re-entrant VT activation patterns.

Conclusion

Our simulation framework has an application in VT RFA intervention planning. It could be used to guide EP explorations and even predict ablation targets pre-operatively. This could reduce intervention duration and improve success rate.

Introduction

Personalized, i.e. patient-specific, simulation of the heart electrical behaviour is a dynamic research field.¹ Building such simulation is often a challenge as mathematical models can be computationally costly and clinical data is sparse and noisy. Most of the work in this area relies on extracting fibrosis from late-enhancement cardiac magnetic resonance imaging (MRI).² In a post-infarct context, this geometry is then classically divided into three different zones based on MRI signal processing: the healthy myocardium, the dense scar, and a grey zone where fibrosis and remaining functional cardiomyocytes coexist. Grouped together, these three zones represent a three-dimensional domain upon which a set of partial differential equations (PDEs) is solved. The model parameters are then chosen differently according to the different zones of the domain to represent the different electrophysiological properties of these tissues. Models typically used are reaction–diffusion models,³ where the reaction part mimics cell-level changes in membrane potential and the diffusion part reflects the ability of cardiac action potential to propagate through specialized structures. Adequately parameterizing such models is a tedious task and solving them takes a lot of computational power and time, incompatible with the classical clinical time frame. These challenges are yet to be overcome for cardiac modelling to enter routine medical practice.

Additionally, before the actual simulation computation, the binary masks (outputs from image segmentation) often need to be transformed into another three-dimensional representation, a three-dimensional mesh, more suited to finite elements PDE solvers.² Choosing adequate parameters for the mesh generation step is critical in order to reach reasonable numerical accuracy and computational time. It requires specific expertise rarely available on day-to-day clinical practice.

Some authors have reported realistic intra-cardiac electrograms simulation,⁴ but the typical outputs of such models are usually visual simulations of the wave front propagation in the cardiac tissue, i.e. activation maps.³ Either way, the goal of such simulations is always to identify abnormal activation patterns linked to fatal arrhythmias. They could improve sudden cardiac death risk prediction,⁵ a major challenge in cardiology.⁶

They can also be used to improve post-infarct re-entrant ventricular tachycardia (VT) radiofrequency ablation (RFA) outcome by contributing to the planning phase of the intervention.^7,–9 In these interventions, the VT is usually induced and mapped (when haemodynamically tolerated) allowing the clinician to identify re-entry isthmus(es) to be ablated. This risky and long targets identification phase still drastically limits the spread and success rate of such interventions. Moreover, most of these patients already carry an implantable cardioverter-defibrillator (ICD), which makes model personalization from MRI difficult if not impossible.

It has been reported that cardiac computed tomography (CT) images could be useful for VT RFA.¹⁰ Computed tomography offers better cost, accessibility, inter-centre and inter-personal reproducibility than MRI. Most notably CTs are not a problem in ICD carriers in whom the image quality and spatial resolution are preserved, as opposed to what happens with MRI. This is important as a majority of scar related VT ablations are conducted in ICD carriers. The better resolution of CT images could also be critical in appreciating scar heterogeneity and thus identifying potential targets.

On CT images, the infarct scar is characterized by an apparent thinning of the myocardial wall,¹⁰ which has been shown to correlate with low voltage areas.¹¹ Some VT re-entry isthmuses are also visible on such images, with a characteristic morphology on myocardial wall thickness maps. It has indeed been recently shown that they appear as a zone of moderate thinning bordered by more intense thinning and joining two areas of normal thickness.¹¹

In this article, we propose a novel framework to simulate activation maps from CT images. It uses CT myocardial wall thickness to parameterize a fast organ-level wave front propagation model in a pipeline visually represented in Figure 1.

Methods

Patient inclusion

The data we used come from seven patients referred for catheter ablation therapy. They all presented chronic (minimum 5 years) infarct and suffered VT episodes despite wearing an ICD.

Cardiac computed tomography images

Images were acquired using a contrast-enhanced electrocardiogram-gated cardiac multi-detector CT using a 64-slice clinical scanner (SOMATOM Definition, Siemens Medical Systems, Forchheim, Germany). Coronary angiography images were acquired during the injection of a 120 mL bolus of iomeprol 400 mg/mL (Bracco, Milan, Italy) at a rate of 4 mL/s. Radiation exposure was typically between 2 and 4 mSv. Images were acquired in supine position, with tube current modulation set to end-diastole. The resulting voxels have a dimension of 0.4 × 0.4 × 1 mm³, superior to what clinically available MRIs can produce. A resulting short-axis slice can be visualized on Figure 1 (left part).

Wall segmentation

On these images, the left ventricle was segmented as follows:

• The endocardial mask was segmented using a region-growing algorithm, the thresholds to discriminate between the blood pool and the wall being optimized from a prior analysis of blood and wall densities.

• The epicardial mask was segmented using a semi-automated tool built within the MUSIC software (IHU Liryc Bordeaux, Inria Sophia Antipolis, France): splines were manually drawn on a few slices and interpolated in-between.

The resulting myocardial wall mask (epicardial mask minus the endocardial mask) has around one million voxels.

Wall thickness estimation

On such images, chronic infarct scars appear as a thinning of the myocardial wall.^10,12 To localize scars, instead of successive dilations of the endocardial masks that were used in other studies,¹¹ we opted for a fully automated and continuous method.¹³ Estimating the thickness of a structure from three-dimensional medical images is not a trivial task. The accuracy can be compromised depending on how the two surfaces are described, and how the distance between them is computed. In this work, we chose to rely on a robust method to compute the shortest path between the endocardium and the epicardium, and then integrate the distance along this path. This method is based on solving a PDE to compute a thickness value at every voxel of the wall mask.¹³ Example results of such thickness maps are presented on Figure 2.

Eikonal model

Numerous models were proposed in the literature to simulate cardiac action potential dynamics. In order to obtain a fast and robust pipeline, we used the Eikonal model of wave front propagation that directly outputs an activation map. It is very robust to changes in image resolution and it can be computed very efficiently with the fast-marching algorithm.

Model inputs are:

• The myocardial wall mask, i.e. the domain of resolution.

• For every voxel of this mask, a local conduction velocity.

• One or more onset point(s), to initiate propagation.

• (Optionally) An artificial unidirectional block, in order to orientate the initial propagation, used to reproduce re-entrant VT patterns.

To solve this equation, we used the fast marching algorithm implementation in the Insight Segmentation and Registration Toolkit (ITK, Kitware, USA). It did not require any meshing step as this implementation directly takes advantage of the regular grid of CT three-dimensional images.

Model parameters

Due to a ‘zig-zag’ course of activation in surviving bundles existing in the infarcted tissue, the depolarization wave front propagates very slowly before finally exiting the scar.¹⁴ We hypothesized that the wave front propagation velocity could be estimated from the wall thickness using a transfer function that would have the following features:

• It should reach a plateau, $vmax$⁠, above physiological wall thicknesses.

• It should be virtually null, $vmin$⁠, below certain thicknesses, where almost no conducting cells remain.

The fast-marching solver needs a ‘stopping value’, i.e. a maximum wave front reaching time. In our experiments, we set it to 500 ms and also assigned this value to non-activated parts of the domain.

In order to model a re-entrant VT pattern, we added an artificial ‘refractory wall’, as previously described by Cedilnik et al.¹⁵ A conduction block orthogonal to the isthmus longitudinal axis was added a few voxels behind the simulation onset point. It allows to model tissue refractoriness, while staying in the simple Eikonal model framework, allowing fast computations and keeping the computational cost very low. This is only needed to reproduce VT activation maps, not for activation maps from controlled pacing.

Automatic ridge detection

In order to detect VT isthmuses based on image criteria alone, we used the objectness filter,¹⁶ implementation of the Insight Toolkit (Kitware, USA). It is a generalization of the Hessian matrix-based vessel enhancement method described by Frangi et al.¹⁷

Briefly, the second order information, i.e. the Hessian Matrix, of an intensity image is used to estimate the probability of each voxel of being part of a tubular structure. In an ideal vessel, the eigenvalues of the Hessian matrix of a voxel are such that:

• The smallest one is close to zero, meaning there is little variation in intensity along a certain direction, along the vessel.

• The other two are significantly negative and almost equal, meaning that in the two direction orthogonal to the vessel direction, there is a brutal change in intensity, i.e. we are going outside of the vessel.

We fed this algorithm with a modified version of the thickness map where only potential isthmuses candidates,¹¹ i.e. moderately thin zones between 2 and 5 mm thickness, remained. We tuned the filter’s parameter to respect the width known to characterize VT isthmuses on CT images.¹¹

We named the output maps obtained this way channelness maps.

This first step leads to a lot of over-detection (false positives) in the scar border zone. To avoid this, we handled the border zone using the following methodology:

• Perform a binary dilation with a radius of 2 mm of the scar zones where the myocardial wall thickness is below 5 mm.

• Perform a binary dilation with a radius of 2 mm of the healthy zones where the myocardial wall thickness is above 5 mm.

• Compute the border zone, defined as the intersection of the two binary masks obtained this way (this intersection is the result of a logical AND operation on them).

• Set the channelness value on the voxel belonging to the border zone to 0.

Comparison with recorded activation maps

In order to gain insight on our model relevance and limitations, we compared our simulations to nine activation maps acquired during RFA interventions in our population of seven patients. The maps were selected when their quality was sufficient to allow proper registration to CT images. These activation maps were produced by the Rhythmia mapping system (Boston Scientific, USA) using a 64-electrodes Orion Catheter. They were recorded during controlled pacing (2 maps) and VT (7 maps). They come in the form of surface triangular meshes in reference frames different from the corresponding CT images.

In order to compare both visually (qualitatively) and numerically (quantitatively) the outputs of our simulations to the reference data, we took advantage of the correlation between low-voltages zones and CT wall thinning areas¹¹ to manually align the two geometries. After fitting an ellipsoid to the point cloud and using the ellipsoid’s centre, we defined each point by spherical coordinates. A visual representation of this registration framework is presented on Figure 4. Finally, it was possible to project the data from the electrophysiology (EP) studies onto the more accurate geometry extracted from the CT images.

Using the ellipsoid long axis and a manually chosen point on the septum, this method also allows the fast and automated generation of bull’s eye plots.

Results

Execution time

The semi-automated segmentation of the left ventricle takes around 10 min for an experienced user of the MUSIC software. The thickness computation takes about 1 min, depending on the number of voxels of the wall mask. The actual simulation takes less than 1 min. The channelness computation is the longest step, around 2–4 min. This informal benchmark results were measured on an Intel i7-5500U-powered laptop.

Automated ridge detection

On all channelness maps, an expert cardiologist and an expert radiologist visually assessed that ridges corresponding to their critical isthmus definition could be detected using the channelness filter (see Figure 5).

Simulated activation maps

In five out of nine maps the general activation patterns of the simulations obtained using our framework match the recordings when using velocities and thicknesses thresholds found in the literature (3/7 VT maps, 2/2 controlled pacing maps). As can be seen on Figure 6, scar late activation is visible both on simulations and on controlled pacing. The characteristic figure of eight activation patterns of re-entrant VT can also be accurately reproduced using our framework, when using a refractory wall to give the propagation an initial direction. While the match between static representations is not always striking, video versions are available as Supplementary material online and offer a better appreciation of dynamics of wave front propagation. This is the usual way recorded activation maps are visualized by the interventional cardiologist in order to decide ablation targets.

For three cases, adapting the thicknesses thresholds could significantly improve the aspect of the simulations (see Figure 7 and discussion below).

Discussion

Image segmentation

In our framework, the simulations rely a lot on how the myocardial wall segmentation is performed. Currently, the endocardial segmentation is robust and perfectly reproducible because it is fully automated. However, the epicardial segmentation is also crucial and is at the moment only semi-automated. Developing a fully automated method for this task is one of our priorities and we expect it to improve our modelling framework’s robustness.

It is worth mentioning that scar transmurality cannot be assessed using the CT scar characterization technique that we used. Using late enhancement CT may be a way to address this issue and we will explore this opportunity in the future.

Model limitations

We voluntarily opted for the simplest cardiac EP model in order to build a fast pipeline. For instance, most of the related work includes an algorithmic fibre orientation estimation that we did not include. We do not believe that it would improve the results a lot, as fibre orientation and its influence are still hard to estimate, particularly in infarcted zones.

In contrast, we believe that an important phenomenon that our model does not address is the front curvature impact. A straight wave front is known to propagate faster than a curved wave front, but this is not the case in our simulations. In some cases, this results in fast, unrealistic, ‘U-turns’ of the excitation front. This could be improved with extensions of the Eikonal model,¹⁸ or by switching to a reaction–diffusion model parameterized on wall thickness. The main challenge will here be to preserve a clinically compatible framework in terms of parameterization complexity and required computational time.

Another limitation of our model is that it only addresses macroscopic re-entrant channel, given the scale at which the simulations are computed. Microscopic re-entries are out of its scope.

Velocities

Our model probably simplifies activation patterns inside the VT isthmuses by globally slowing down the propagation in them. A recent animal study¹⁹ indeed suggests that the depolarization slows down when entering and exiting the isthmus but propagates normally in between. The wave front would then be slowed down inside the core of the isthmus. However, we think that our global slowing down in the whole isthmus is a relevant averaging simplification at the organ scale. In the case where our further work would show the pro-eminence of this phenomenon, an option would be to use channelness maps to accordingly change velocities in isthmuses.

Thickness thresholds

Based on data from a recent study,¹¹ we defined global thickness thresholds across all patients to discriminate between healthy and scar zones. Specifically, we set them to below 2 mm for ‘pure scar’ and above 5 mm for healthy myocardium. However, in some VT cases, the wave front propagation crossed scar zones too fast and matching the recording proved more challenging.

We noticed that in those cases, for instance the one shown in Figure 7, the simulation output could be drastically improved by shifting these thresholds to higher values, i.e. by raising the $p$ parameter of our thickness to velocity transfer function. Using patient-specific thresholds both for thickness maps analysis and model parameterization could improve the use of CT images in the context of VT RFA. These patient-specific thresholds could possibly be determined by the disease history or image features such as the distribution of thickness for a specific myocardial wall.

Parameter optimization

Model

Using the registration of image and EP data in the common frame of reference described in the Methods section, it was possible to define a metric to quantify differences between simulations and ground truths. Using the mean square difference between these data, we therefore optimized the thickness to velocity transfer function parameters (⁠$vmax$⁠, $vmin$⁠, $p,$ and $r$⁠, see Figure 3) to minimize the discrepancy between simulations and recordings. To this end, we used the covariance matrix-evolution strategy algorithm,²⁰ suited for this type of non-convex optimization problem. While minimizing the mean square difference did improve the numerical matching of the simulations to the recording, the overall activation pattern actually was often degraded. This is partly due to noise in the target EP data and to the fact that the mean square difference metric may not be adapted to preserving the activation patterns like the typical figure of eight shape of re-entrant VT. How to define this pattern-aware metric is still an open question.

Channelness

The output of the channelness is also highly depending on setting the right parameters for the objectness mapping, including a channelness threshold. These parameters need to be optimized and validated on a larger multi-centric database.

Conclusion: a clinically compatible modelling framework

The simulation framework presented in this article may be useful for EP modelling to enter clinical practice. It requires little human interaction and expertise to go from imaging data to EP simulations. It is very fast and would even be faster by developing implementations more suited to such purpose than the generic ones we currently use. The availability and reproducibility offered by CT images and such automation could help personalized cardiac modelling to enter clinical applications.

Funding

The research leading to these results has received European funding from the ERC starting grant ECSTATIC (715093), the Equipex MUSIC (ANR-11-EQPX-0030), and the IHU LIRYC (ANR-10- IAHU-04).

Conflict of interest: none declared.

References