A rule-based method for predicting the electrical activation of the heart with cardiac resynchronization therapy from non-invasive clinical data.

BACKGROUND: Cardiac Resynchronization Therapy (CRT) is one of the few effective treatments for heart failure patients with ventricular dyssynchrony. The pacing location of the left ventricle is indicated as a determinant of CRT outcome.
OBJECTIVE: Patient specific computational models allow the activation pattern following CRT implant to be predicted and this may be used to optimize CRT lead placement.
METHODS: In this study, the effects of heterogeneous cardiac substrate (scar, fast endocardial conduction, slow septal conduction, functional block) on accurately predicting the electrical activation of the LV epicardium were tested to determine the minimal detail required to create a rule based model of cardiac electrophysiology. Non-invasive clinical data (CT or CMR images and 12 lead ECG) from eighteen patients from two centers were used to investigate the models.
RESULTS: Validation with invasive electro-anatomical mapping data identified that computer models with fast endocardial conduction were able to predict the electrical activation with a mean distance errors of 9.2 ± 0.5 mm (CMR data) or (CT data) 7.5 ± 0.7 mm.
CONCLUSION: This study identified a simple rule-based fast endocardial conduction model, built using non-invasive clinical data that can be used to rapidly and robustly predict the electrical activation of the heart. Pre-procedural prediction of the latest electrically activating region to identify the optimal LV pacing site could potentially be a useful clinical planning tool for CRT procedures.

Introduction

Cardiac resynchronization therapy (CRT) has emerged as an effective therapy for heart failure patients with ventricular conduction disturbances, such as left bundle branch block (LBBB), resulting in a dyssynchronous ventricular activation. CRT aims to restore this dyssynchronous activation and the location of the left ventricular (LV) pacing lead has repeatedly been identified as an important contributor to patient response (Singh et al., 2011; Ypenburg et al., 2008).

The optimal pacing site has been proposed to be located at the latest activating regions (Khan et al., 2012; Saba et al., 2013; Zanon et al., 2014) and outside of scar (Khan et al., 2009; Leyva et al., 2011; Singh et al., 2011). However, this approach has two short comings. Firstly, in cases where scar tissue causes a circumferential block in the activation wave, the latest region to activate can be located adjacent to scar. Secondly, the optimal location is based on the state of the patient prior to implant and is determined without consideration of the location of the right ventricle (RV) pacing lead, which may alter the location of the latest activating region (Mafi Rad et al., 2014). Recent studies have found that increased rate of endocardial activation patterns (Sohal et al., 2015) and reduced QRS duration (Derval et al., 2014) result in an improved response to CRT, this could potentially provide a better cost function for optimizing optimal pacing location.

As the heart is paced from the RV and LV electrodes that can be positioned independently, arguably pacing locations need to be optimized not on the intrinsic activation pattern in the heart prior to CRT, but taking into account the effect of both leads on the activation pattern once CRT is implanted. Optimizing therapy planning using the activation pattern post implant requires the ability to predict how different pacing configurations will alter the activation pattern for given lead locations.

Computer models of cardiac electrophysiology have been developed to perform patient specific simulations to predict electrical activation patterns. This approach has been used for predicting changes in activation patterns in CRT and to investigate the optimal atrioventricular delay, interventricular delay and lead locations (Lee et al., 2017; Miri et al., 2009; Reumann et al., 2007a, 2007b; Villongco et al., 2016).

The muscle fiber directions within the heart are important in determining the pattern of electrical activation throughout the ventricles, as the electrical activation spreads predominantly along the direction of the myofibers. Measurements of the fiber angles (Holmes et al., 2000; Scollan et al., 1998; Seemann et al., 2006) have found the fiber orientation to be relatively consistent in the ventricles and well described by rule-based methods (Bayer et al., 2012). Therefore, in recent years, studies have adopted rule-based fiber methods to approximate the fiber directions in electrophysiology models of human hearts (Aguado-Sierra et al., 2011; Crozier et al., 2015; Kayvanpour et al., 2015; Kerckhoffs et al., 2012; Lee et al., 2017; Okada et al., 2017; Sermesant et al., 2012; Tobon-Gomez et al., 2013; Villongco et al., 2016).

Previous models have personalized ventricular activation and predicted activation patterns for CRT using invasive measures of the electrical properties of the heart, such as electro-anatomical mapping (EAM) of the epicardium via the coronary veins or the endocardium by contact or non-contact mapping (Aguado-Sierra et al., 2011; Crozier et al., 2015; Lee et al., 2017; Sermesant et al., 2012; Tobon-Gomez et al., 2013; Villongco et al., 2016). These studies required invasive data, collected intra-procedurally, meaning that these models are unable to be used for prospective prediction of the electrical activation of the heart or to track the response to CRT.

In recent studies, non-invasive ECG data from conventional 12 lead ECGs or from body surface potential mapping (with commercial versions recording up to 256 ECGs) were used to parameterize electrophysiology models (Giffard-Roisin et al., 2017; Kayvanpour et al., 2015; Okada et al., 2017). In these studies, the cardiac substrate was assumed to be homogeneous with uniform conduction velocity (Giffard-Roisin et al., 2017) or with fast endocardial conductivity (Kayvanpour et al., 2015; Okada et al., 2017). These studies used non-invasive data to constrain the models, however these modeling frameworks required complex and iterative fitting to parameterize the onset site/s of the electrical activation (Giffard-Roisin et al., 2017; Okada et al., 2017) and the conductivities of the myocardium (Giffard-Roisin et al., 2017) and the fast conducting endocardium (Kayvanpour et al., 2015; Okada et al., 2017).

A rule-based method for defining the electrophysiology parameters would be useful in facilitating the development of electrophysiology models in individual patients. However, the material properties of the heart have been less well described, with experimental studies reporting a range of conduction velocities (0.07–0.75 m/s) along the fiber direction, with a 2–7 fold increase in electrical conduction velocity along longitudinal direction in comparison to the transverse direction of the myofibers (Clerc, 1976; Roberts et al., 1979; Roth, 1997).

Experimental studies have also found heterogeneous conduction velocities throughout the ventricles. In addition to the three to six fold increase in the conduction velocity in the Purkinje system (Draper and Mya-Tu, 1959; Myerburg et al., 1978), the sub-endocardial layer has also been found to have a 1.5–4 fold increased electrical conduction velocity in comparison to the normal myocardium (Draper and Mya-Tu, 1959; Myerburg et al., 1978; Strik et al., 2011). In experimental studies of the canine heart, the spread of the electrical activation in the endocardium was found to start simultaneously across the RV and starting in the bottom third of the LV and spreading in an apical to basal direction in the LV (Myerburg et al., 1972), while in human hearts activation in the LV began from endocardial areas in the lower third of the free wall, mid septal wall and in the basal area underneath the mitral valve (Durrer et al., 1970). A number of electrophysiology simulation studies have thus incorporated a fast endocardial layer representing the Purkinje network or the increased conduction velocity in the sub-endocardial layer (Aguado-Sierra et al., 2011; Hyde et al., 2015; Kayvanpour et al., 2015; Okada et al., 2017).

Additional heterogeneities have been reported to influence activation. The conduction velocity of scarred myocardium have been found to be reduced in comparison to normal tissue (de Bakker et al., 1993; Peters et al., 1993). Clinical studies have also found that the electrical activation pattern in the LV can be classed as type I/II in patients suffering from LBBB, with non-contact mapping studies showing a C-shaped or U-shaped pattern of electrical activation where there is a region of functional block in the anterior or posterior regions of the LV (Auricchio et al., 2004; Jia et al., 2006). A decrease in the conduction velocity has also been observed in the septum in canines and humans with LBBB (Prinzen and Auricchio, 2008; Sohal et al., 2014; Strik et al., 2013). These heterogeneities have been reported across different pathologies and species, however their relative importance in predicting the electrical activation pattern has not been determined.

In this study we investigated the impact of scar, slow septal conduction, fast endocardial conduction and the presence of type I/II activation patterns on predicting cardiac electrophysiology for CRT patients using simple models of tissue heterogeneity. The simplest rule-based methods for defining the electrical heterogeneities present in the ventricles needed to robustly predict the activation patterns based on non-invasive imaging of the patient anatomy and ECG measurements was identified.

Methods

Study population

This study was conducted on data collected from fourteen patients who underwent a de novo CRT device implantation at the Maastricht University Medical Center and underwent a cardiac magnetic resonance (CMR) imaging before implantation. An additional four RV-pacing patients upgrading to CRT were recruited at Guys and St Thomas’ Hospital. All patients had a I or II CRT indication according to the ESC guidelines (LV ejection fractions of < 35% with mild to severe heart failure symptoms with LBBB according to specific criteria or non-LBBB with a QRS duration of > 150 ms). Data was collected from patients at the two centers with the study protocols approved by the respective local ethics committees.

Non-invasive data

In the CMR data set, images were acquired with a 1.5–3T Philips scanner (Achieva/Ingenia/Intera) prior to device implantation. Steady state free precession CMR sequence (slice thickness 6–8 mm, field of view 320–384 mm, matrix 256–560 × 256–560) was used to image the heart in 2-chamber, 3-chamber, 4-chamber views as well as a short axis stack covering the entire LV. Contrast-enhanced (CE) CMR images were also acquired with a 2D gradient echo inversion recovery sequence (slice thickness 7–10 mm, field of view 300–384 mm, matrix 240–576 × 240–576). The CMR images were manually segmented using customized software (CAAS MRV3.4, Pie Medical Imaging) to personalize the cardiac anatomy. Areas of cardiac infarct were semi-automatically segmented from the CE CMR images with the full-width half maximum method (Nguyen et al., 2017).

In the four CT cases, 3D whole heart CT images (slice thickness 0.5 mm, field of view 185–280 mm x 185–280 mm x 118–191 mm, and in-plane resolution 0.36–0.54 mm) were acquired with a Philips iCT256 or Siemens Somatom Force scanner. RV paced patients were scanned and subsequently underwent a CRT upgrade procedure, where a LV lead is implanted into the coronary sinus (CS) for epicardial pacing. The Philips model-based automatic segmentation tool (Peters et al., 2007) was used to segment the ventricles and the main CS branch. Semi-automatic segmentation of the coronary branches was then performed using CRT segmentation software (Mountney et al., 2017).

Non-invasive measures of the electrical activation of the heart were obtained with 12 lead ECG during RV pacing. The QRS onset and duration from one of the frontal leads were used to estimate the start of the electrical activation and the total time of activation in the ventricles.

Electroanatomical mapping

At both centers, 3D EAM of the coronary veins was performed using Ensite NavX (St Jude Medical) intra-procedurally (Mafi Rad et al., 2014). In the de novo patient cases, an RV lead was implanted in the RV apex region, guided with fluoroscopy X-ray images and identified on the EAM. A unipolar sensing and pacing guidewire connected to the Ensite NaxX system was inserted into the CS and extended into the CS tributaries, to generate a 3D map of the CS and its tributaries as well as measuring the local activation times (LAT) of these sites with RV pacing. Across the 14 cases a large epicardial area of the LV free wall was mapped, as shown in Appendix A. The LAT, defined as the duration from onset of the QRS complex on the surface ECG to the steepest downslope on the intracardiac unipolar electrogram, was recorded throughout the procedure and was calculated as a percentage of the QRS duration. In the CMR data set, the EAM 3D anatomy was mapped to the segmentations of the CS using the fusion function of the Ensite NavX system as a post processing step by a clinician using the X-ray fluoroscopy images obtained during the procedure as guidance.

Electrical activation models

CMR cases

Patient specific electrophysiology models were created from the 14 CMR data sets. The anatomy for each case was personalized from the segmentations. The location of the RV lead was identified from the EAM mapped onto the models using the fusion function in the Ensite NavX system. The locations where the LAT was measured with the guidewire were projected onto the closest point on the LV epicardial surface of the patient specific meshes with a root mean square error of 3.0–9.1 mm for the 14 CMR cases.

A rule-based method was used to determine the fiber orientations (Bayer et al., 2012). The propagation of the electrical activation from pacing at the RV apex was simulated using an Reaction-Eikonal model (Neic et al., 2017) using the Cardiac Arrhythmia Research Package (CARP) (Niederer et al., 2011; Vigmond et al., 2003). The Eikonal conductivity term was fitted to match the QRS data using a grid search of low cost eikonal activation simulations. The Reaction-Eikonal conductivity was then input into the Reaction-Eikonal equations. This results in a predicted activation time error of no more than 1%. The conduction velocity was characterized by assuming the QRS duration was equal to the total activation time across the ventricles, with the computed QRS duration set to be within a tolerance of 5 ms. The conduction velocity in the ventricles was assumed to be transversely anisotropic, with the conduction velocity transverse to the fiber directions set to be 40% of the conduction velocity along the myofiber directions (Sermesant et al., 2012). The importance of including scar, slow septal conduction, fast endocardial conduction and the presence of type I/II activation patterns in predicting epicardial activation patterns were also explored using additional electrophysiology models [Fig. 1].

Fig. 1

Six electrophysiology models were investigated with the inclusion of scar, functional block in the anterior or posterior regions of the LV, slow septal conduction, and fast endocardial conduction used to determine the importance of these factors in accurately simulating the electrical propagation across the ventricles.

In the scar regions and functional block regions, the electrical activation was assumed to be completely blocked. While the slow septal conduction velocity was assumed to be 50% of the normal myocardium (Strik et al., 2013). A 1 mm thick fast endocardium layer was assumed to have a six-fold increased conduction velocity along the fiber direction compared to the normal myocardium (Draper and Mya-Tu, 1959; Myerburg et al., 1978).

CT cases

Four patients were receiving an upgrade from RV pacing to CRT and were not eligible for CMR. CT whole heart images were segmented using Philips Model Based Segmentation Framework and used to create patient specific models. The location of the RV pacing lead was identified on the CT images. The coronary venous EAM was recorded intra-procedurally using a unipolar sensing and pacing guidewire connected to the Ensite NavX system. The anatomical 3D map was registered to the segmented coronary sinus from the CT images using CMISS (www.cmiss.org) host mesh fitting, with clinical guidance based on X-ray fluoroscopy imaging to determine which branches were mapped. As the segmentations of the CS and the patient specific meshes were derived from the same CT image, the mapped coronary venous EAM were also mapped onto the meshes. The locations where the electrical delay was measured with the guidewire were projected onto the closest point on the LV epicardial surface of the patient specific meshes with a root mean square error of 3.1–8.6 mm for the 4 CT cases.

The LAT from RV pacing was clinically measured with EAM in the coronary sinus tributaries. The models were evaluated by comparing the simulated LAT with intra-procedural clinical LAT measurements. The relative LAT error was taken as the difference between the model simulations and the clinical measurements. The temporal error for each simulation was calculated by multiplying the relative LAT error by the QRS duration for each case. Positive temporal error values indicate that the model activates slower than the clinical measurements, while negative values indicate that the model activates faster. The distance error measurement for each simulation was calculated as the minimum distance between the location of the EAM measurement projected onto the myocardium and the nearest tissue with an activation time that matches the measured activation time.

Results

The electrical activation for each of the models was simulated for each of the 14 CMR cases with the conduction velocity along the fiber directions for the bulk myocardium characterized by the QRS duration. The propagation of the electrical activity from pacing at the RV for the six simple electrophysiology models (normal, with scar, with anterior or posterior block, with slow septal conduction and with fast endocardial conduction) were simulated as shown in Fig. 2. The LAT errors were computed for each model in each case as shown in Fig. 3.

Fig. 2

RV pacing with a six-fold increase in the ventricular endocardium was simulated for the 14CMR cases, with the electrical activation time normalized as a percentage of the QRS duration. The latest activated site is highlighted with a white dot.

Fig. 3

The simulated local activation time normalized as a percentage of the QRS duration (LAT) was compared against the clinical LAT for each of the measured sites in the coronary sinus venous branches for the model with fast endocardial conduction.

The resulting temporal and distance errors for the six electrophysiology models are shown in Fig. 4. The results and analysis in the next section are from a leave-one-out cross validation, and the reported scores are the average of the mean temporal and distance errors. One-way ANOVA was used to compare the temporal and distance errors between the six simple electrophysiology models, and it was found that there were statistically significant differences in the temporal and distance errors for the models (p-value < 0.001). Tukey post-hoc tests indicated that only 2 models had significantly different temporal means (posterior functional block model (–5.1 ± 1.2 ms) and the fast endocardial conduction model (–7.0 ± 1.2 ms)), compared to the rest of the models (Normal: 7.6 ± 1.2 ms, Scar: 7.4 ± 1.2 ms, slow septum: 8.0 ± 1.2 ms, anterior functional block: 6.6 ± 1.2 ms), however the absolute temporal error mean values were comparable. Tukey post-hoc test indicated that the fast endocardial conduction model was the only model that had a significantly reduced mean distance error (9.2 mm ± 0.5 mm) in comparison to the other models (15.6–16.9 mm ± 0.5 mm) (see Appendix C).

Fig. 4

Simulations were run for the 14 CMR cases, with a basic model, inclusion of scar (for cases 9–14), slow septum, fast endocardial conduction, anterior or posterior functional block. Boxplots of the (a) temporal error, (b) distance error, and (c) conduction velocities are shown for each model. The conduction velocities which are between the physiologically plausible range (0.07–0.75 m/s) are highlighted in green.

To further constrain the model selection, several criteria were imposed on the solution of the models to ensure the simulations reflect known physiological constraints. In experimental studies, the conduction velocity of the ventricular myocardium has been found to range between 0.07–0.75 m/s (Clerc, 1976; Draper and Mya-Tu, 1959; Kleber and Rudy, 2004; Roberts et al., 1979; Sano et al., 1959; Weidmann, 1970), the models which extended outside this range were excluded from consideration. Endocardial mapping in LBBB cases have found the latest point of activation to lie close to the lateral wall of the LV (Auricchio et al., 2004), thus we also excluded models where the latest site of activation was in the septum. The comparison of model conduction velocity and physiological constraints is shown in Fig. 4C. The only model which fulfilled both criteria was the model with fast endocardial conduction. The mean temporal error for the models with six-fold fast endocardial conduction in the CMR cases was −7.0 ± 0.3 ms, mean distance error was 9.2 mm ± 0.1 mm, with mean model conduction velocity of 0.33 ± 0.08 m/s.

Sensitivity to model parameters

The sensitivity of the distance error measure to changes in the anatomy, fibre orientation, anisotropy ratio, and the slow septal conductivity were also investigated. It was found that regardless of changes in the anatomy, fibre orientation, anisotropy ratio, slow septal conductivity, the conclusion that fast endocardial conduction was the most important factor remained consistent (see Appendix B for details).

In addition to models with a six-fold increased conduction velocity in the endocardium, simulations were also run for up to ten-fold (Aguado-Sierra et al., 2011) increases in the endocardial conduction velocity (see Appendix B3). It was found that as the fast endocardial conduction ratio increased from 1.0 (normal model), the distance error between the model simulations and the clinical measurements gradually improved. We observed that there were no significant differences in the mean distance errors for five-fold to ten-fold increases in the FEC ratio (7.9–10.0 mm ± 0.4 mm). While the myocardial conduction velocity for all the fast endocardial conduction models remained within the values reported in experimental studies.

Experimental canine and human studies have shown that the normal depolarization of the heart proceeds in an apical to basal direction, with the greatest concentration of arborizations of the Purkinje network in the lower third of the heart in the apical-basal direction (Myerburg et al., 1972; Spach et al., 1963). Additional simulations were carried out with setting the fast endocardial conduction to be in the lower third of the heart. One-way ANOVA was used to compare the models with varying fast endocardial conduction ratios and it was found that, when the fast endocardial conduction was only located on the lower third of the heart, there were no differences in the mean distance error measures compared to the normal model (see Appendix B3).

Cross-modality validation

A set of 4 CT patient cases was used for cross-modality validation of the rule based material property models (normal, with anterior or posterior block, with slow septal conduction and with fast endocardial conduction) (Fig. 5). No scar data was segmented, so the scar model was not tested with these CT data sets.

Fig. 5

(a) CT images were acquired for 4 CRT upgrade patients and were automatically segmented using the Philips model based tool. (b) The EAM from Ensite NavX (white) was mapped to the coronary sinus semi-automatically segmented from the CT images (black). (c) The electrical activation across the ventricles from RV pacing was simulated. (d) The EAM local activation times were mapped onto the model (spheres) and the RMS distance error of the mapped EAM sites onto the models is stated for each case. (e) The LAT as a percentage of the QRS duration was compared between the clinical measurements and the model simulations with a six-fold increased conduction velocity in the ventricular endocardium.

The temporal and distance error measures were calculated for the six electrophysiology models (Fig. 6a and b). One-way ANOVA indicated that there were no significant differences in the mean temporal errors (p-value > 0.05), while there were found to be significant differences in the mean distance errors of the six electrophysiology models (p-value < 0.001). Tukey post-hoc tests indicated that six-fold fast endocardial conduction (7.5 ± 0.7 mm) had a significantly improved mean distance error measure compared to the other models (Normal: 31.4 ± 0.7 mm; with anterior (32.2 ± 0.7 mm) or posterior block (27.5 ± 0.7 mm), and with slow septal conduction (32.0 ± 0.7 mm). The model fitted conduction velocities were within physiologically plausible ranges (0.07–0.75 m/s) for the 4 CT data sets (0.43 ± 0.18 m/s) only for the fast endocardial conduction models (Fig. 6c).

Fig. 6

Simulations were run for the 4 CT cases, with a basic model, slow septum, fast endocardial conduction, anterior or posterior functional block. Boxplots of the (a) temporal error, (b) distance error, and (c) conduction velocities are shown for each model. The conduction velocities which are between the physiologically plausible range (0.07–0.75 m/s) are highlighted in green.

Discussion

In this study we developed a model that robustly simulates the electrical activation across the heart using non-invasive data. Rule-based methods for defining the electrical properties of the heart were investigated for 14 CMR cases and an additional 4 CT patient cases. The importance of the inclusion of scar, functional block, slow septal conduction and fast endocardial conduction in simulating the electrical activity across the heart was studied.

Our results indicate that fast endocardial conduction is the most important factor in constraining the model to physiologically plausible parameters consistent with experimental studies (latest electrical activation in a non-septal region with conduction velocities of 0.07–0.75 m/s) and resulting in sub-10 mm errors on average. This mean distance error needs to be considered in the clinical context. For CRT pacing lead optimization, the activation patterns needs to identify the optimal AHA segment. Endocardially and epicardially each AHA segment in the eighteen models had a mean size of 31.0 ± 4.6 mm x 14.4 ± 1.4 mm and 45.2 ± 4.3 mm x 30.2 ± 4.0 mm, respectively. The mean distance errors of the CMR cases (9.2 ± 0.5 mm) and the CT cases (7.5 ± 0.7 mm) are well within the sensitivity to identify the optimal AHA segment.

The underlying mechanism behind fast endocardial conduction is still unknown, with either retrograde activation of the Purkinje network or intrinsic properties of the endocardial layer hypothesized as reasons for an increase in the endocardial conductivity (Akar et al., 2004; Berruezo et al., 2004; Hyde et al., 2015; Strik et al., 2012; Taccardi et al., 2008). Experimental studies have found a transmural gradient in cellular properties across the ventricles, with increases in the cell areas (Stoker et al., 1982), gap junction density (Glukhov et al., 2012), and sodium channel densities (Gaborit et al., 2007) corresponding to increased conduction velocity in the endocardium in comparison to the epicardium. However, these experimentally observed changes in the endocardial cells are unable to fully explain the reported 1.5–4-fold increase in the conduction velocity (Draper and Mya-Tu, 1959; Myerburg et al., 1978; Strik et al., 2011).

An alternative explanation for fast endocardial conduction is retrograde activation of the Purkinje network, which has been observed during RV pacing (Damato et al., 1970). In experimental studies, Purkinje fibers arborizations have primarily been found in the lower third of the ventricles (Myerburg et al., 1972; Spach et al., 1963), though endocardial breakthrough sites have also been identified further up near the base of the left ventricle (Durrer et al., 1970). Simulations were run for fast endocardial conduction in the RV and the lower third of the LV, however it was found that the mean distance error had no significant improvements in comparison to the normal model, regardless of the increase in endocardial conduction velocity, unlike the models with a fast endocardial layer extending fully up the ventricles where the mean distant error improved as the fast endocardial conduction ration increased (Appendix B).

To date, there is a lack of evidence to clearly differentiate between the two mechanisms for fast endocardial conduction. In animal studies, the Purkinje network has been found to have a three to six fold increase in the conduction velocity in comparison to the myocardial tissue (Draper and Mya-Tu, 1959; Myerburg et al., 1978; Sano et al., 1959), while the superficial layers of the sub-endocardial ventricular muscle was found to have a 1.5–4 fold increase (Draper and Mya-Tu, 1959; Myerburg et al., 1978; Strik et al., 2011). Our simulations were also unable to distiguish between the two mechanisms, with five-fold to ten-fold increases in the fast endocardial conduction yielding similar temporal and distance errors for models where the fast endocardial conduction layer was defined as extending from the apex to base of the ventricles.

Other simulation studies of cardiac electrophysiology have specifically included in a Purkinje network represented as a 1D cable network (Romero et al., 2010; Vigmond and Stuyvers, 2016) or as activation sites (Kerckhoffs et al., 2012; Tobon-Gomez et al., 2013) based on experimental studies (Durrer et al., 1970; Usyk et al., 2002). The Purkinje system anatomy is individual specific, however in-vivo Purkinje fibers are unable to be imaged with current methods. Electrophysiology simulation studies have used manually delineated models of the Purkinje network (Romero et al., 2010) or randomly generating activation onset sites on the endocardial surface (Tobon-Gomez et al., 2013), while patient-specific modeling of the Purkinje network has relied on complex iterative methods and invasive electro-anatomical data (Cárdenes et al., 2015; Vergara et al., 2014). The aim of this study was to investigate simple rule-based methods that can be parameterized from non-invasive clinical measurements, so we have assumed that a fast endocardial layer approximates either the Purkinje network or the rapid conduction of the sub-endocardial layer. Further investigation into the importance of Purkinje network activation in comparison to this approximation remains to be done.

In contrast to clinical studies where the scar plays a significant role in response to CRT, our study has shown that scar was not the main determinant of activation times. However, the conclusions in this study need to be taken in the context of the aims of this paper, which were to investigate the relative importance of various factors effecting the electrical substrate of the ventricles. Scar, in addition to its effects on the electrical properties, has also been found to effect the mechanical properties of the cardiac tissues, with a higher density of collagen fibers in scarred tissue leading to an increase in the mechanical stiffness and decreased contractility (Jugdutt et al., 1996). Although in this study, we have found that scar is not important, it is in the context of the electrophysiology pattern, while the mechanical response of the heart is not modeled. In addition, in the study simulations a major effect of the scar would only be seen if it is completely transmural, otherwise the simulated electrical activation would simply go around the scar. So in a clinical context, it would be important to know where the scar is, but in the simulation of the electrical activation pattern, it is less important than including the fast endocardial conduction.

Limitations

The electrical activation in the heart has been shown to spread primarily along the fiber directions, so the myofiber orientations have a significant effect on the electrophysiology in the heart (Clerc, 1976; Roberts et al., 1979; Roth, 1997). In recent years, there has been increasing interest in personalizing the fiber orientation of the heart with in-vivo (Dou et al., 2003; Gamper et al., 2007; Geerts et al., 2002; Nguyen et al., 2014; Nielles-Vallespin et al., 2013; Toussaint et al., 2010, 2013; Wei et al., 2015) or ex-vivo (Helm et al., 2005; Holmes et al., 2000; Rohmer et al., 2007) diffusion tensor MRI (DTMRI). Though DTMRI is a promising technology for personalizing the fiber orientations in the heart, current challenges regarding the signal-to-noise ratio, long image acquisition times, low resolution and motion artifacts need to be overcome prior to using this method to personalize the fiber orientations in electrophysiology models of the heart. In previous computer modeling studies on rat (Bishop et al., 2009) and canine (Bayer et al., 2012) models, DTMRI-based and rule-based methods for determining the fiber orientation yielded similar electrical activation patterns. In this manuscript, we have thus assumed that the fiber orientations can be described adequately with rule-based methods.

Activation time measurements were only recorded on the free wall of the LV epicardium from within the coronary veins. This places inherent limitations on the locations that activation times can be measured. To confirm that measurements were made from across the epicardium, we have mapped recording sites onto an AHA plot (see Appendix A) and found multiple measurements points from all LV epicardium regions. The model was therefore only evaluated against epicardial activation time measurements. While complete activation time mapping would be ideal, patient measurements must be clinically tractable. We focused our measurements on the epicardium as this is the location of conventional CRT pacing and we would like to predict the latest epicardial activation time, using simulations, to identify potential optimal CRT pacing lead locations.

Another limitation of the models is that the regions of the heart were categorized into normal bulk myocardium, scar, functional block, slow septal conduction and fast endocardial conduction layers, with the electrical modeling parameters assumed to be homogeneous within each region. The conduction velocity in each region was set to be zero (scar or functional block) or scaled to be relative to the bulk myocardial tissue. However, regional heterogeneities have been found even across normal cardiac tissue, with varying properties across the heart transmurally, apical-basally and between the ventricles. In addition the infarcted regions were segmented based on MRI images, with no consideration of the heterogeneities across the border zone of the scar data. These additional heterogeneities could affect the conduction velocities across the heart. However, the degree of heterogeneity in the cardiac substrate is still unknown, particularly in diseased hearts. The myocardial conduction velocity was parameterized based on the QRS duration, a singular measure from the non-invasive 12 lead ECG traces. Other information can be derived from the ECG traces, such as the electrical axis which describes the main direction of the electrical depolarization in the heart. This information could potentially be used to further personalize the non-bulk myocardial conductivities in the heart model, such as differentiation of the RV and LV endocardial conduction velocities (Kayvanpour et al., 2015). However, the accuracy and robustness of these methods would still need to be investigated in future studies.

Clinical implications

In this study, we have developed and validated a novel approach to generate a realistic endocardial and epicardial electroanatomic activation map using only non-invasive modalities which are commonly used in standard clinical care (CMR/CT and 12-lead ECG). This approach was shown to work across 2 imaging modalities (CMR/CT). As only non-invasive data is required to predict the electrical activation of the ventricles, these simulations can be run prior to CRT implantation. In contrast to previous work relying on bidomain or monodomain models, our simulations used Eikonal equations to approximate the spread of the electrical activation across the heart, which requires a fraction of the computational costs (Neic et al., 2017). In this study, each simulation of the heart electrical activation took < 30 s on a standard desktop computer. While body surface potential mapping is a commercially available non-invasive technique that can be also be used to map the electrical activation of the heart, its cost remains a barrier to widespread clinical adoption. Our method presented here is fast, cheap and uses non-invasive clinical images and data acquired as part of standard clinical practice allowing for easy integration into standard clinical workflows to prospectively predict the latest point of electrical activation in the heart prior to device implantation for CRT patients.

Conclusion

The combination of a rule-based method for the fiber orientation (Bayer et al., 2012) with a rule-based description of the heterogeneous material properties of the electrical system in the heart allows for electrophysiology models to be developed with minimal non-invasive data (medical images such as CT or MRI and 12 lead ECG to determine the QRS duration) prior to device implantation. There is evidence that the optimal LV pacing site for CRT lies in the latest electrically activating region (Zanon et al., 2014). A model that incorporates a fast endocardial conduction layer is able predict the electrical activation to sub-10 mm accuracy. These model can potentially be a useful clinical planning tool to prospectively predict the latest point of electrical activation in the heart prior to device implantation (Table A9).