Accuracy of prediction of infarct-related arrhythmic circuits from image-based models reconstructed from low and high resolution MRI

Identification of optimal ablation sites in hearts with infarct-related ventricular tachycardia (VT) remains difficult to achieve with the current catheter-based mapping techniques. Limitations arise from the ambiguities in determining the reentrant pathways location(s). The goal of this study was to develop experimentally validated, individualized computer models of infarcted swine hearts, reconstructed from high-resolution ex-vivo MRI and to examine the accuracy of the reentrant circuit location prediction when models of the same hearts are instead reconstructed from low clinical-resolution MRI scans. To achieve this goal, we utilized retrospective data obtained from four pigs ~10 weeks post infarction that underwent VT induction via programmed stimulation and epicardial activation mapping via a multielectrode epicardial sock. After the experiment, high-resolution ex-vivo MRI with late gadolinium enhancement was acquired. The Hi-res images were downsampled into two lower resolutions (Med-res and Low-res) in order to replicate image quality obtainable in the clinic. The images were segmented and models were reconstructed from the three image stacks for each pig heart. VT induction similar to what was performed in the experiment was simulated. Results of the reconstructions showed that the geometry of the ventricles including the infarct could be accurately obtained from Med-res and Low-res images. Simulation results demonstrated that induced VTs in the Med-res and Low-res models were located close to those in Hi-res models. Importantly, all models, regardless of image resolution, accurately predicted the VT morphology and circuit location induced in the experiment. These results demonstrate that MRI-based computer models of hearts with ischemic cardiomyopathy could provide a unique opportunity to predict and analyze VT resulting for from specific infarct architecture, and thus may assist in clinical decisions to identify and ablate the reentrant circuit(s).


INTRODUCTION
Myocardial infarction (MI), a condition characterized by reduced viability of cardiac myocardium due to insufficient blood supply, is a leading cause of lethal ventricular tachyarrhythmia (VT) worldwide (Go et al., 2014). A promising therapeutic approach is to identify and eradicate the arrhythmogenic substrate through an invasive procedure known as catheter ablation (Delacretaz and Stevenson, 2001). This technique utilizes catheters in an electrophysiological (EP) study that records surface electrical activity to identify the location(s) of the VT reentrant circuits, which are then targeted for ablation. However, an EP study is limited to interrogation of the cardiac surfaces alone, which fails to take into account the complex 3D architecture of the infarct and could lead to missing intramural arrhythmogenic substrates. Additionally, the point by point nature of the technique results in long procedure times that increase risk for complications in patients (Stevenson et al., 2008). These limitations have resulted in a success rate of only 58% for first time VT ablation procedures (Callans et al., 1998). Clearly, there remains a need for novel methodologies to safely and effectively identify the specific MI regions that harbor VT reentrant pathways which can then guide successful post-MI VT ablation.
Recent advances in non-invasive imaging techniques have shown promise in supplementing the information gained from EP studies by providing a priori information on infarct structure (Ashikaga et al., 2007;Dickfeld et al., 2008;Tian et al., 2010;Rutherford et al., 2012). In particular, contrast-enhanced (CE) MRI has become the gold standard for assessing the location, transmurality, and composition of MI (Ashikaga et al., 2007;Perez-David et al., 2011;Fernández-Armenta et al., 2013). Recently, the presence of peri-infarct or border zone surrounding the necrotic scar, also known as gray zone (GZ) based on their appearance as regions of intermediate intensity in the CE-MRI scans, have been shown to correlate with increased risk of post-MI mortality (Yan et al., 2006), increased likelihood of spontaneous clinical VT (Roes et al., 2009), and increased VT inducibility via programmed stimulation (Schmidt et al., 2007). Although CE-MRI is a powerful tool to visualize the complexity of the MI structure, it does not provide any insight into the electrical activity in the heart, particularly the location of the VT reentrant circuits. Additionally, while targeting GZ for ablation could result in increased procedure success, the complex distribution of multiple GZs throughout the heart makes identification of the specific arrhythmogenic substrate difficult (Nayyar et al., 2014).
The limitations of EP studies and non-invasive imaging can be addressed by recent advances in individualized, computational modeling of the heart. Recently, our group and others have demonstrated that models of post-MI hearts reconstructed from high resolution ex-vivo MRI can be used to gain insight into the mechanisms of post-MI VT initiation and maintenance (Pop et al., 2011;Rantner et al., 2012;Arevalo et al., 2013). Importantly, these studies confirmed that specific regions of GZ within the complex infarct structure promote reentry initiation and maintenance. Work has also already begun on bringing the technology closer to clinical applicability by reconstructing models from low resolution clinical MRI (Ng et al., 2012;Ashikaga et al., 2013;Ringenberg et al., 2014). However, difficulty in obtaining electrical recordings during VT has made it difficult to fully assess the accuracy of computational models in predicting the morphology and location(s) of VT reentrant circuits. Additionally, it is not known whether models generated from low resolution MRI, with its loss of infarct geometry detail, can still accurately predict post-MI VT.
Thus, the goal of this study is (1) to develop experimentally validated, individualized computer models of infarcted pig hearts based on high resolution images that can accurately predict post-infarction VT morphology and reentrant circuit location and (2) to examine the accuracy of the reentrant circuit location prediction when models of the same hearts are instead reconstructed from low clinical-resolution MRI scans. Determining this accuracy will provide the needed information for the development of patient-specific models from clinical MRI scans and their usilization in predictive simulations regarding infarct-related arrhythmia treatment.

Electrophysiological Data and Image Acquisition
Part of the research in this study utilized experimental data and CE-MRI previously published in the study by Ashikaga et al. (2007). In that study, six domestic pigs underwent MI induction via occlusion of the left anterior descending artery. Approximately 10 weeks post-MI, bipolar electrograms were recorded from 300 to 380 ventricular epicardial sites via a multielectrode sock. VT was then induced via programmed stimulation delivered to a bipolar electrode randomly placed on the epicardium. After the procedure, the hearts were explanted and high resolution ex-vivo CE-MRI was acquired. The final image matrix was 256 × 256 × 256, and the image resolution was 0.39 × 0.39 × 0.39 mm 3 . In two of the six pig hearts, the induced VT had chaotic and heterogeneous patterns. These two pig data were excluded from the present study.
In order to mimic typical clinical MRI resolutions, original images were downsampled by eliminating the high-frequency components in the Fourier space of the original images using a 3D rectangular mask. The size of the rectangular mask was calculated as a function of the desired resolution in each dimension (Nishimura, 2010). Following this, an inverse FFT was applied to reconstruct the downsampled images. As a result, the corresponding original and the downsampled images had the same voxel (and matrix) size for the segmentation process. Two levels of downsampling were performed resulting in medium resolution (1.5 × 1.5 × 4 mm 3 ) and low resolution (1.5 × 1.5 × 8 mm 3 ) image stacks.

Hi-res Images
The Otsu threshold method, available in the software package Seg3D (SCI Institute, Salt Lake City, UT, USA), was used to automatically segment the myocardium from the surrounding bath (Figure 1, top row). The myocardium was then further segmented into three tissue types, normal tissue, gray zone (GZ), and scar, using a modified gray-level thresholding algorithm (Ng et al., 2012;Arevalo et al., 2013). First, the region of infarct was manually selected within each slice. Next, the mean and standard deviation (SD) of the voxel intensities in the region outside the infarct were calculated. The gray-level value corresponding to mean + 4 SD was assigned as the upper limit for normal tissue (normal max ). All voxels with intensities less than normal max were classified as normal tissue. Finally, the range of intensities within the infarct tissue (infarct range ) was calculated by obtaining the difference between the maximum intensity within the infarct region (infarct max ) and normal max . Voxels with intensity between normal max and (normal max + 50% of infarct range ) were classified as GZ. Voxels with intensity more than normal max + 50% of infarct range were classified as scar.

Med-res and Low-res Images
The downsampled image stacks were segmented using a previously described segmentation and interpolation method (Prakosa et al., 2014;Ukwatta et al., 2015). Briefly, CardioViz3D (INRIA, Sophia-Antipolis, France), was used to seed points and obtain the surfaces for the epi-and endocardial boundaries of the ventricles (Figure 1, bottom row). The same gray-level thresholding method described in the hi-res image processing section was used to further segment the myocardium. Due to the overall lower range of pixel intensities resulting from the downsampling, the gray-level value corresponding to normal max was calculated from mean + 3 SD as opposed to mean + 4 SD. Calculation of the gray-level thresholds for GZ and scar was done as previously described.

Mesh Generation and Fiber Assignment
The finite element computational meshes were constructed directly from the segmented images using a previously described approach (Prassl et al., 2009). The procedure preserved the fine geometric details of the ventricles and the different infarct zones. Finally, fiber orientations were assigned in the mesh using a rule-based method (Bayer et al., 2012).

Modeling of Electrophysiological Properties
Mathematical description of cardiac tissue was based on the monodomain representation (Plank et al., 2008). The scar was modeled as passive tissue. Anisotropic conductivities in the normal myocardium were assigned values of 0.28 S/m in the longitudinal myofiber direction and 0.026 S/m in the transverse direction, which resulted in conduction velocity value of 45 cm/s in the longitudinal direction, matching swine conduction velocities reported by Vetter et al. (2005). Within the GZ, the transverse conductivity was decreased by 90% to represent connexin 43 remodeling, as has been done previously .
Membrane behavior in the normal tissue was represented by the Luo-Rudy dynamic model (LRd) (Luo and Rudy, 1994). The same membrane model was used, with modifications based on experimental data, to represent the electrophysiology of GZ cells, as we have previously described , namely: peak sodium current was decreased to 38% of the original value in the LRD model, peak L-type calcium current was decreased to 31% of the original value, peak potassium currents I Kr and I Ks were decreased to 30 and 20% of the original values, respectively. The modifications resulted in a decreased action potential upstroke velocity, decreased peak action potential magnitude, and increased action potential duration. Previous research by our team  has indicated that representation of the GZ as a homogeneously remodeled region provides an adequate approximation in the calculation of the morphology and location of the infarct-related reentrant circuits.

Simulation Protocol and Data Analysis
All simulations were performed using the software package CARP (CardioSolv, LLC) on a parallel computing platform (Vigmond et al., 2003;Plank et al., 2008). To examine the arrhythmogenic propensity of the post-MI ventricular models, programmed electrical stimulation (PES), similar to the protocol used in the pig experimental study, was simulated (Ashikaga et al., 2007). The models were paced from an epicardial location for 6 beats (S1) at a cycle length of 600 ms followed by a premature stimulus (S2) initially given at 90% of S1 cycle length. The timing between S1 and S2 was progressively shortened until VT was induced. If VT was not induced, a second premature stimulus (S3) was delivered after S2. If VT was not induced either, a third premature stimulus (S4) was delivered after S3. PES was performed at 10-12 epicardial sites to fully assess VT inducibility. Three seconds of VT was simulated.
Pseudo-ECGs were calculated in order to facilitate comparisons of global activation patterns during VT (Bishop and Plank, 2011). Two points separated by 20 cm were placed near the base of the heart with the line connecting them perpendicular to the base-apex axis. The difference in the extracellular potential was then calculated between the two points. To quantify difference between two ECG waveforms, the mean absolute deviation (MAD) metric was computed.
Where X and Y are the pseudo-ECG waveforms of length n, and X and Y are the average values of the waveform. MAD varies between 0 and 100%, corresponding to identical and completely different waveforms, respectively.
Finally, the difference in location of the reentrant circuits induced in the different models was quantified. VT reentrant circuit location was estimated from the activation maps. In the case of reentrant circuits located on the epicardial or endocardial surface, the midpoint of the scar that anchors the reentry was identified (Figures 2-5, yellow star). If the epicardial VT pattern had a breakthrough morphology, intramural slices were used to locate the reentrant circuit and identify the scar that anchors the circuit.  Table 1 summarizes the volumes of the normal myocardium, GZ, and scar in the reconstructed hearts. A paired Student's t-test revealed that there was no significance difference between total ventricular volumes between the Hi-res models (154.4 ± 7.5 cm 3 ) and either of the Med-res or Low-res models (155.3 ± 9.4 and 157.9 ± 11.5 cm 3 , respectively). However, in all four hearts, the models reconstructed from the Med-res and Low-res images had significantly larger GZ and scar volumes compared to the Hi-res model (p < 0.05). Looking at the composition of each model, the Med-res and Low-res models had a significantly smaller percentage of normal myocardium (90.8 ± 0.8% and 92.0 ± 1.1%) compared to the Hi-res model (93.5 ± 1.5%). This demonstrated that the increase in GZ and scar volumes was due to expansion of the infarct onto normal tissue. This was to be expected since downsampling results in partial volume effect at the infarct periphery.

Accuracy of VT Simulations
For Pig 1, the ventricular geometry reconstructed from the different resolution images is presented in Figure 2A, 1st row. The overall structure of the anterior infarct in the Med-res and Low-res models matched well the Hi-res model. However, detailed structures in the infarct periphery, particularly near the apex, were lost as the image resolution decreased.
Despite the loss in infarct structure detail, the induced VT activation maps remained similar among the three models (Figure 2A, 2nd row). In the Hi-res model, PES resulted in the induction of a VT that manifested as breakthrough activity on the epicardium at the RV-septum junction of the anterior wall. The same epicardial breakthrough activity was induced in the Med-res and Low-res models. In the Hi-res model, a transmural slice through the epicardial breakthrough site (dashed line) revealed that the reentrant pathway driving VT was anchored to intramural scar located in the antero-septal region (Figure 2A, 3rd row). In the Med-res and Low-res models, the reentrant circuit was anchored to scar located on the LV-side of the septum endocardium. Despite the difference in the reentrant pathways, the distance between the scar that anchored the reentrant circuit in the Hi-res model and the scars that anchored the reentrant circuit in the Med-res and Lowres models was relatively close (13.6 and 9.4 mm, respectively, Table 2).
The pseudo-ECGs, further demonstrated the accuracy of the downsampled models in predicting post-MI VT (Figure 2B). The representative traces during 1 VT cycle showed that the cycle length and overall morphology recorded from the three models were similar. This was evident in the low MAD score when the Hi-res trace was compared with the Med-res and Low-res traces (11.8% and 11.0%, respectively).
Finally, all three models accurately predicted the experimentally induced post-MI VT morphology. In the experiment, Pig 1 was inducible for VT with epicardial breakthrough (Figure 2C). The breakthrough was located on the anterior wall at the RV-septum border toward the base, similar to what was induced in all three models. The simulation results for Pig 2 followed the same trend as the Pig 1 results. The 1st row of Figure 3A shows that the reconstructed models from the downsampled images accurately resolved the ventricular and infarct geometry. In all models, simulated PES resulted in induction of VT with breakthrough activity on the anterior epicardium ( Figure 3A, 2nd row). In both the Hi-res and Med-res models, the underlying reentrant activity was an intramural circuit anchored to scar located in the septum (Figure 3A, 3rd row). The reentry-anchoring scar in the Med-res model was located 14.5 mm from the scar that anchored the reentry in the Hi-res model ( Table 2). The close match in the reentrant patterns also resulted in a small MAD score (10.5%, Figure 3B). For the Low-res model, the breakthrough activity was driven by a reentrant circuit on the RV septum endocardial wall. The organizing center of this reentrant circuit was located 18.3 mm from the reentry-anchoring scar in the Hires model. The difference in activation pattern compared to the Hi-res model resulted in a slightly higher MAD score (17.8%) compared to the Med-res model. All three models successfully predicted the epicardial VT activation pattern and location of breakthrough site that was experimentally recorded in the same pig ( Figure 3C). Figure 4; the ventricular and infarct geometry of the Hi-res model were again successfully resolved in the models reconstructed from the downsampled images ( Figure 4A, 1st row). In the Hi-res model, PES resulted in the induction of breakthrough activity on the epicardium at the RV-septum junction (Figure 4A, 2nd row). The transmural view reveals that the underlying reentrant circuit was confined within the septum at the RV apex ( Figure 4A, 2nd row). The epicardial activation map of the VT induced in the Med-res model, had a breakthrough pattern with the site of first epicardial activation identical to the one induced in the Hi-res model. The reentrant driver for this VT was located on the RV septum; whose organizing center was located 6.2 mm from the reentrant driver in the Hi-res model ( Table 2). The similarity in the morphology  and location of the reentrant circuit resulted in a low MAD score of 11.3% ( Figure 4B). Meanwhile, in the Low-res model, PES resulted in the formation of a figure-of-eight reentrant circuit on the epicardium that manifested as breakthrough on the RV endocardium. While the VT morphology differed to what was observed in the Hi-res model, the Low-res model's reentrant circuit was still closely located to the reentrant circuit induced in the Hi-res model (7.9 mm). Interestingly, the simulated VT in the Pig 3 Hi-res model did not match the experimentally induced VT.

Results for Pig 3 are presented in
In the experiments, Pig 3 was inducible for VT with a figure-of-8 epicardial reentry morphology ( Figure 4C); the experimental pattern was matched best by the lLow-res model. Regardless, the location of the experimental reentrant circuit was close to the epicardial breakthrough and reentry site simulated in the Hi-res and Med-res models. Due to the difference in activation patterns, the Low-res model had a higher MAD score (21.1%) when compared to the Hi-res model. Finally, the results for Pig 4 continue the trend of accurate geometric match between the three models observed in the other pigs ( Figure 5A, 1st row). In this pig, PES in the Hi-res model resulted in the induction of a figure-of-eight reentrant pattern on the epicardium located on the anterior RV-septal junction toward the apex (Figure 5A, 2nd row) and manifests as breakthrough activity on the LV endocardial surface (Figure 5A, 3rd row). On the other hand, the Med-res and Low-res models resulted in the induction of breakthrough activity on the epicardium with the epicardial breakthrough located near the figure-of-8 reentry site induced in the Hi-res model. In both Med-res and Lowres models, the epicardial breakthrough was driven by rotors located on the LV endocardium. The Med-res and Low-res rotors were located closely to the High-res rotor (13.0 and 9.4 mm, respectively). Due to difference in the VT morphologies between the Hi-res model and the Med-res and Low-res models, the MAD score for Pig 4 was higher than the ones calculated for the other pigs (28.8 and 32.0%, respectively, Figure 5A, 3rd row).
The experimental recording revealed that Pig 4 was inducible for figure-of-8 reentry located at the anterior apex ( Figure 5C). The simulation results of the Hi-res model successfully matched the experimentally induced VT pattern and location. Meanwhile the Med-res and Low-res models were unsuccessful in accurately predicting the experimental VT morphology, but the two models were still able to accurately identify the location of the arrhythmogenic substrate in the heart.

DISCUSSION
In this work, we successfully created image-based models of post-MI ventricles reconstructed from both high-resolution exvivo MRI and as well as from clinical MRI scans. Simulations with these image-based models were used to predict the post-MI arrhythmogenic reentrant circuits following PES. While there have been previous attempts in simulating post-MI electrophysiological activity (Vigmond et al., 2009;Pop et al., 2011;Ng et al., 2012), this is the first study to use experimentally obtained VT activation maps to validate model predictions as well as demonstrate that the locations of the in-silico predicted VT circuits are fairly well preserved in models reconstructed from low resolution MRI. These results pave the way for the translation of in-silico modeling of post-MI arrhythmogenesis into clinical utility.
Previous experimental and clinical studies have identified the heterogeneous structure of infarcted tissue as the arrhythmogenic substrate in post-MI hearts. Specifically, studies have shown that the presence of partially viable tissue interdigitated with electrically inert scar promote conduction slowing and provide pathways that could harbor reentrant circuits (Stein et al., 2008;Suk et al., 2008). The experiments that produced the retrospective pig data used in the current study were the first to combine high resolution ex-vivo CE-MRI with multi-electrode sock recording in order to determine the relationship between VT circuits and specific infarct structure (Ashikaga et al., 2007). The results demonstrated that islands of viable myocardium within the scar correlated with epicardial breakthrough sites in the case of centrifugal VT pattern and with the isthmus in the case of epicardial reentry VT pattern. While this study hinted at the utility of MRI-based scar analysis to predict the location of the arrhythmogenic substrate, the image resolution for exvivo preparations is currently not achievable in the clinical setting.
We demonstrate in this study that the complexity of infarct structure can be reasonably reconstructed from low resolution images. In all four pigs, the models reconstructed from both the Med-res and Low-res images had infarct geometries similar to that obtained from the Hi-res ex-vivo images. In the infarct periphery, downsampling of the images led to partial volume effect leading slight overestimation of the infarct volume as compared to the Hi-res images. There was also difficulty in resolving infarct structure close to the apex. This was due to uneven diffusion of gadolinium toward the apex resulting in decreased maximum image intensity close to the apex (Schelbert et al., 2010). Despite this loss in infarct structural detail, the overall infarct architecture was preserved in the downsampled images. Additionally, the percent composition of scar and GZ in the infarct was comparable between the Hi-res images and the downsampled images.
Simulations with the reconstructed hearts demonstrated the accuracy of the models in predicting the location of infarctrelated reentrant circuits. In the Hi-res models, three of the four hearts successfully predicted the epicardial morphology of the experimentally obtained VT in the ex-vivo hearts. In all four hearts, the simulations successfully predicted the location of the reentrant circuit or the epicardial breakthrough site. In the three hearts with epicardial breakthrough morphology, the simulations revealed that the underlying reentrant circuit was located intramurally or on the endocardium. This was consistent with the Ashikaga et al's findings that swine hearts were more likely to have endocardial than epicardial reentrant circuits (Ashikaga et al., 2007).
Using the simulation predictions from the Hi-res models as ground truth, the predicted VT location remained fairly accurate for the models reconstructed from the downsampled images.
The VT activation maps in the Med-res and Low-res models matched the morphology of the epicardial activation recorded in 75 and 50%, respectively, of the Hi-res pig hearts. In order to quantitatively measure the difference of electrical conduction between different images resolution, we use MAD score of the pseudo-ECGs (Vadakkumpadan et al., 2012). The higher overall accuracy of the Med-res models compared to the Low-res models was reflected in the lower mean MAD scores (15.6 ± 7.6 vs. 20.5 ± 7.6%). In any case, the relatively low MAD scores for both Med-res and Low-res models demonstrate the accuracy of the downsampled models in simulating overall VT activation patterns.
The downsampled models had a higher success in predicting the location of the reentrant circuits. The predicted reentrant circuits in the Med-res and Low-res models were closely located from the circuits induced in the Hi-res models (11.8 ± 3.3 vs. 11.3 ± 4.1 mm). Since current radiofrequency ablation lesions can have up to 14 mm width (Petersen et al., 1999), the predicted circuits in the downsampled models were within one ablation lesion away from the "true" VT reentrant circuit location predicted in the Hi-res models. Since clinical ablation usually involve the delivery of multiple lesions in one area, using models reconstructed from clinical images to guide VT ablation appears feasible.
The results of this study could have implications for the potential use of patient-specific simulations of infarct-related VT in predicting the optimal sites of VT ablation. The relatively high level of preservation of the locations of the organizing centers of reentry or sites of breakthrough upon downsampling of the images, as found here, indicates that it might be feasible to use patient-specific heart models reconstructed from clinical MR scans to determine the VT morphology in the given patient and thus to provide guidance as to what ablation lesions would best terminate VT. Future studies are needed to address this subject of major clinical significance.

LIMITATIONS
A limitation of the current study is the small sample size of only four pig hearts. The original experimental study included more pigs but a few had chaotic epicardial VT patterns, indicative of complex intramural reentrant activity, which was difficult to match in silico. Several of the pigs underwent endocardial basket electrode recording which was difficult to register with the MRI scan. Purkinje fibers were also not included in the models which might have affected the propagation patterns during VT, especially at sites distal from the reentry origin. Finally, some of the ex-vivo images had signal noise due to gadolinium leak. In these cases, the images had to be manually cleaned.