Competing Risks in Patients with Primary Prevention Implantable Cardioverters Defibrillators: Global Electrical Heterogeneity and Clinical Outcomes (GEHCO) Study

Study design

The study design and protocol has been previously described.¹⁰ In brief, we conducted a retrospective, multicenter cohort study that included data from six academic medical centers in the United States: OHSU in Portland, OR, Veteran Administration Portland Healthcare System (Portland VA) in Portland, OR, Beth Israel Deaconess Medical Center (BIDMC) in Boston, MA, the University of Colorado (U Colorado) in Aurora, CO, Cedars-Sinai Medical Center (CSMC) in Los Angeles, CA, and Stanford University (Stanford U) in Palo Alto, CA. We included all patients above the age of 18 years with chronic systolic heart failure (HF) due to infarct-related CM and/or NICM, who underwent implantation of any ICD [including single chamber, dual-chamber, cardiac resynchronization therapy defibrillator (CRT-D), or subcutaneous] for primary prevention of SCD between the years of 1996 and 2019, and had an available digital 12-lead ECG recording. Patients with inherited channelopathies and cardiomyopathies, congenital heart disease, and those with ICD implanted for secondary prevention were excluded.

ECG analysis: global electrical heterogeneity measurement

A digital 12-lead ECG recorded around the time of initial ICD implantation was collected (median 0 days (interquartile range 0-29 days). Raw, digital ECG signals were analyzed in the Tereshchenko laboratory at OHSU, as previously described.^{3, 7, 11, 12} At least two physician-investigators manually labeled each cardiac beat (SH, EB, LGT). The Kors transformation matrix was used to transform 12-lead ECGs into orthogonal XYZ vectorcardiograms.¹³ The time-coherent global median beat was constructed using the dominant type of beat present in the ECG, and the origin point of the vectorcardiogram was identified.¹² We included three categories of median beats: normal (N) median beats included normal sinus rhythm, atrial pacing, junctional rhythm, and ectopic atrial rhythm; ventricular paced (VP) median beats included right ventricular or biventricular pacing; supraventricular (S) median beats included atrial fibrillation or atrial flutter with a consistent QRS morphology. QT interval was corrected using Bazett, Fridericia, and Hodges equations. The presence of premature ventricular complexes (PVCs) on 12-lead ECG was noted.

We constructed spatial peak and area QRS and T vectors as well as their vector sum (SVG) and measured their magnitude, direction (azimuth (orientation in the XZ plane) and elevation (orientation in the XY plane)), and spatial QRS-T angles (Figure 1).^{7, 11, 12} Scalar approximations of the SVG were measured via the sum absolute QRST integral (SAIQRST)^14-16 and the QT integral on vector magnitude (VMQTi).¹¹ Quality control of automated ECG analysis was performed to verify appropriate beat identification and fiducial point annotation (KTH). Open-source MATLAB (MathWorks, Natick, MA, USA) code for GEH measurement is provided at https://physionet.org/physiotools/geh & https://github.com/Tereshchenkolab/Origin.

• Download figure
• Open in new tab

Figure 1.

GEG measurements. A. Transformation of a 12-lead ECG to orthogonal XYZ ECG. B. The spatial ventricular gradient (SVG) is a vector formed from the sum of the area QRS and area T wave vectors. The QRST angle is the 3-dimensional angle between the area QRS and area T wave vectors. C. Orientation of angles in 3D space. Positive angles are posterior and negative angles are anterior. D. SAIQRST is the sum of the absolute areas of the X, Y, and Z leads. E. The vector magnitude lead (VM) is calculated as VM=√(X²+Y²+Z²). The VM QT integral (VMQTi) is the area under the QT ECG signal on VM lead.

Primary Competing Outcomes

The primary outcome was defined as the first sustained ventricular tachyarrhythmia event (VT or VF) with appropriate ICD therapy (either anti-tachycardia pacing or shock). All-cause death without preceding appropriate ICD therapy served as the primary competing outcome.

Clinical characteristics and ICD programming

Patients’ demographics, clinical history, and laboratory values at the pre-implant clinical evaluation were abstracted from the medical record. Estimated glomerular filtration rate (eGFR) was calculated using the Chronic Kidney Disease Epidemiology Collaboration equation.¹⁷ ICD detection parameters were obtained at the time of device implant or at the time of appropriate ICD therapies if they occurred during follow-up. To harmonize VT/VF detection duration programming across different device manufacturers, we grouped patients¹⁸ in (1) long detection duration (number of intervals to detect (NID) VT ≥ 30/40 or VT detection duration ≥ 10s or VF detection duration ≥ 3s) or (2) short detection duration (NID VT <30/40 or VT detection duration < 10s or VF detection duration < 3s) categories (Supplemental Table 1).

Statistical analyses

Normally distributed continuous variables were reported as means and standard deviation (SD). Circular variables were reported as mean and 95% confidence interval (CI). An unadjusted comparison of circular variables (spatial QRS-T angle, SVG azimuth, and SVG elevation) was performed using the Wheeler-Watson-Mardia test. ANOVA (for normally distributed continuous variables) and Pearson’s χ² test (for categorical variables) were used to compare ECG characteristics in patients with and without primary competing outcomes. Unadjusted analyses were performed using the dataset with missing values.

As we conducted a retrospective study, the data contained missing values for reasons unrelated to the study itself. Missingness data for individual covariates are reported in Supplemental Table 2. The most likely reason for missing values (incomplete data in the medical record) was unrelated to the unobserved data but likely was dependent upon the observed data (e.g., ICD therapy), which is characteristic for data missing at random. Therefore, we employed multiple imputations using chained equations (MICE)¹⁹, which does not require the multivariate normal assumption. Our study’s investigation of data missingness showed an arbitrary pattern. Observations with missing primary outcomes were excluded. We used predictive mean matching with 20 nearest neighbors for continuous variables, logistic regression for binary variables, ordinal logistic regression for ordinal variables, and multinomial logistic regression for nominal variables with 66 imputations, exceeding the minimum required by the von Hippel method.²⁰ Convergence was confirmed. Relative variance increase, fraction missing information, and relative efficiency were assessed for each multivariable model and confirmed acceptable.

After examining a normal quantile plot, we included QRS-T and SVG elevation angles in all adjusted statistical analyses without transformation because their distributions were nearly normal. The SVG azimuth angle is ranging from -180° to +180°. As recommended,^21-23 we transformed SVG azimuth by doubling its value and then adding 360°.

As the risk of sustained VT/VF leading to appropriate ICD therapy competes with the risk of death without appropriate ICD therapy, we employed cause-specific hazards functions, estimated using Cox proportional hazards models. For a model with sustained VT/VF and appropriate ICD therapy outcome, all-cause death event was censored at the date of death if there was no documented sustained ventricular tachyarrhythmia at any time since implant until death. Accordingly, for a model with death without preceding appropriate ICD therapy outcome, sustained VT/VF with appropriate ICD therapy event was censored at the arrhythmic event date. Alive and event-free patients were censored at the date of the last device follow-up.

We compared the strength of association of demographic, clinical, ECG, and device characteristics between the two competing outcomes (appropriate ICD therapies and all-cause death) by testing the null hypothesis that the coefficients for a given variable (b₁ and b₂) were the same across the two competing outcomes by calculating Z-scores: where s.e. is a standard error, and determining their statistical significance.

To study an association of GEH with the primary competing outcomes, we constructed four models. All models (1-4) with GEH exposure variables were adjusted for the type of median beat, mean RR interval, presence of PVCs, and possible distortion of median beat by premature beats. Model 1 was adjusted for demographic characteristics (age, sex, race) and study center. Model 2 included model 1 variables with additional adjustment for known clinical risk factors of the primary competing outcomes [LVEF, NYHA class, CM type, hypertension, diabetes, stroke, atrial fibrillation, eGFR_CKD-EPI, and medications (use of beta-blockers, ACEi/ARBs, and class 1 or III antiarrhythmic drugs (AAD)]. Model 3 contained all covariates included in model 2, and, in addition, device characteristics: device type (ICD or CRT-D), manufacturer, and programming (VT and VF zone cutoffs and detection duration category). Model 4, in addition to model 3 variables, was adjusted for traditional ECG metrics: heart rate, QRS duration, and Hodges-corrected QT interval. To standardize comparisons, we expressed continuous variables per 1 SD.

To investigate previously reported discrepancy in the association of SAIQRST with ventricular tachyarrhythmias/SCD between studies that included entirely or mostly infarct-related CM (ICM) patients^{3, 24} and studies that included mixed ICM and NICM patients,^{16, 25} we tested the interaction of CM type with SAIQRST and VMQTi.

To appropriately consider nonlinear (e.g., U-shaped) hazards, associations of continuous GEH variables with primary competing outcomes were also studied using adjusted (model 1) Cox regression models incorporating cubic splines with 4 knots.

To test the study findings’ robustness under the missing at random assumption, we conducted a sensitivity analysis using an imputed dataset based on complete ECG data. The dataset without missing pre-implant ECG data included 2251 patients; 480 of them had sustained VT/VF with appropriate ICD therapy, and 401 died without appropriate ICD therapy.

Statistical analyses were performed using STATA MP 16.1 (StataCorp LP, College Station, TX). A P value < 0.05 was considered statistically significant. STATA do files are available at https://github.com/Tereshchenkolab/gehco.

Study Population

After excluding the ineligible patient records and those with missing outcomes, our study population included 2,668 patients (Figure 2). Clinical characteristics are shown in Table 1.

• Download figure
• Open in new tab

Figure 2.

Study flowchart.

Over a median retrospective follow-up of 4 years, 541 patients experienced sustained ventricular tachyarrhythmia leading to appropriate ICD therapy [incidence rate 50.1 (95% CI 46.0-54.5) per 1000 person-years], and 479 patients died without appropriate ICD therapy [incidence rate 44.0 (95%CI 40.2-48.1) per 1000 person-years]. The average rate of appropriate ICD therapy was similar to the rate of death without appropriate ICD therapy: ∼ 5% per year.

An average heart rate, QRS duration, and QTc were within the normal range (Table 2). Patients who died without appropriate ICD therapy had a higher heart rate than those who remained alive and event-free (P=0.030). Patients who received appropriate ICD therapies were more likely to have PVC recorded on 12-lead ECG, smaller SVG magnitude and SAIQRST, wider QRS-T angle, and SVG vector pointing more superiorly, in comparison to events-free patients.

Demographic and clinical characteristics associated with primary competing outcomes

In minimally adjusted model 1 (Supplemental Table 3), male sex, black race, atrial fibrillation, use of AADs, history of MI and revascularization procedures, NYHA class III, and higher BUN were associated with increased risk of appropriate ICD therapies. In contrast, greater age, eGFR, and LVEF were associated with a lower risk of appropriate ICD therapies. Adjustment for cardiac disease substrate, risk factors (model 2), and device programming (model 3) explained the association of race, atrial fibrillation, MI and revascularization, NYHA class, and BUN with sustained VT/VF, which became non-significant. In fully adjusted analysis (model 4), only age, male sex, use of AADs, LVEF, and eGFR remained significantly associated with sustained VT/VF and appropriate ICD therapies (Supplemental Table 3).

In model 1 assessing the competing outcome of death without appropriate ICD therapy (Supplemental Table 4), age, ICM, MI, revascularization history, NYHA class II-IV, diabetes, atrial fibrillation, and BUN were associated with a higher risk of all-cause death, while use of beta-blockers, ACEIs, higher LVEF, and higher eGFR were associated with a lower risk of death. Adjustment for the HF characteristics and risk factors (model 2), device characteristics (model 3), and traditional ECG metrics (model 4) only slightly attenuated an association of clinical risk factors with the competing mortality outcome.

Even after comprehensive adjustment (model 3), several clinical risk factors had significantly different cause-specific hazards associated with the two competing outcomes (Supplemental Table 5). For every 13 years (1 SD) increase in age, the hazard of appropriate ICD therapy decreased by 20% (Figure 3), whereas the hazard of death without appropriate ICD therapy increased by 32%. Diabetes increased the hazard of death without appropriate ICD therapy by 45% but was not associated with appropriate ICD therapy. As expected, NYHA class and BUN had stronger associations with death than with appropriate ICD therapy. The use of beta-blockers and ACEIs was associated with a reduced hazard of death but not with the hazard of appropriate ICD therapy.

• Download figure
• Open in new tab

Figure 3.

Statistically significantly different (P<0.05) or similar (P≥0.05) adjusted (model 3) cause-specific Cox hazard ratios and 95% confidence interval of two competing outcomes: sustained ventricular tachyarrhythmia with appropriate ICD therapies (red ovals) and all-cause death without prior ventricular tachyarrhythmia with appropriate ICD therapies (blue rectangles).

In contrast, there was no statistically significant difference in adjusted cause-specific hazards of two competing outcomes for male sex, ischemic CM, atrial fibrillation, the use of AADs, LVEF, and eGFR (Figure 3).

Device characteristics associated with primary competing outcomes

In minimally adjusted model 1, ICD device type and programming were significantly associated with sustained VT/VF and appropriate ICD therapies (Supplemental Table 3). Adjustment for HF characteristics, risk factors, comorbidities, and traditional ECG metrics (models 2-4) did not change the strength of the association. In contrast, very few device characteristics were associated with mortality without appropriate ICD therapy (Supplemental Table 4).

In adjusted competing risk analysis, device characteristics had significantly different cause-specific hazards associated with the two competing outcomes (Supplemental Table 5). Unsurprisingly, for each additional 15 bpm increase in VT zone threshold, the hazard of sustained ventricular tachyarrhythmias with appropriate ICD therapies decreased by 25% (Figure 3), but it did not affect all-cause mortality without sustained VT/VF. Consistent with prior observations, more aggressive ICD programming was associated with a 56% higher risk of appropriate ICD therapies but was not associated with the competing death outcome.

CRT-D implantation was associated with a lower risk of both competing outcomes (Figure 3).

Association of GEH with primary competing outcomes

In minimally adjusted analysis, the presence of PVCs and all GEH metrics were associated with sustained VT/VF and appropriate ICD therapies (Supplemental Table 3). However, the adjustment for known clinical risk factors and device characteristics explained that association for SAIQRST and VMQTi. Notably, after full adjustment in model 4 (including heart rate, QTc, and QRS duration), area and peak spatial QRS-T angle and SVG elevation, area SVG magnitude, and peak SVG azimuth were associated with appropriate ICD therapies.

In comparison, only resting heart rate, spatial QRS-T angle, and area SVG magnitude associated with the competing mortality outcome. The strength of associations diminished from model 1 to model 4, suggesting confounding by the HF substrate (Supplemental Table 4).

In minimally adjusted competing risk analysis, heart rate and QTc had significantly different cause-specific hazards associated with the two competing outcomes (Supplemental Table 5). However, these differences were explained by differences in clinical risk factors. In adjusted models, QRS duration and corrected QT interval were not associated with either appropriate ICD therapies or total mortality. After adjustment, only PVCs and peak SVG azimuth had significantly different cause-specific hazards associated with two competing outcomes (Figure 3). The presence of PVCs was associated with a 42% higher risk of appropriate ICD therapies but did not associate with competing mortality risk.

Each additional 1 SD increase in peak SVG azimuth (+36°, backward) was associated with a 12% increase in the hazard of appropriate ICD therapies and a 10% reduction in the hazard of total mortality without appropriate ICD therapy. Accordingly, 1 SD reduction in peak SVG azimuth (−36°, forward) was associated with a 12% decrease in the hazard of appropriate ICD therapies and a 10% increase in the risk of death without appropriate ICD therapy. Figure 4 illustrates opposing directions of hazards of two competing outcomes across peak SVG azimuth distribution. Area SVG azimuth showed a U-shaped association with appropriate ICD therapy (Supplemental Table 3) but not competing death risk (Supplemental Table 4).

• Download figure
• Open in new tab

Figure 4.

Adjusted (model 1) risk of (A) sustained VT/VF with appropriate ICD therapy and (B) competing death without appropriate ICD therapy associated with peak SVG azimuth. Restricted cubic spline with 95% confidence interval shows a change in cause-specific Cox hazard ratio (HR) (Y-axis) in response to peak SVG azimuth change (X-axis). The 50^th percentile of the SVG azimuth is selected as a reference. Knots of the peak SVG azimuth are at (−21)-42-75-137 degrees.

Spatial QRS-T angle and SVG magnitude were associated with both appropriate ICD therapy and death without appropriate ICD therapy. For each 1 SD (33°) increase in area QRS-T angle, we observed an ∼20% increase in the hazard of both outcomes (Figure 3). Each additional 1 SD increase in area SVG magnitude (25 mV*ms) was associated with a 25% reduction in the hazard of appropriate ICD therapies and a 14% reduction in the hazard of death without appropriate ICD therapy. For area SVG elevation, each 1 SD (31°) increase was associated with a 14% increase in the hazard of appropriate ICD therapies.

Sensitivity analysis (Supplemental Table 6) showed similar results.

Interaction of SAIQRST with cardiomyopathy type

Cardiomyopathy type significantly modified the association of SAIQRST (and VMQTi) with the primary arrhythmic (but not competing death) outcome (Supplemental Table 7). After adjustment for HF characteristics, risk factors, and comorbidities (model 3), for each additional 111 mVms increase in SAIQRST, the hazard of sustained VT/VF with appropriate ICD therapies increased by 29% in ICM patients, but decreased by 22% in NICM subgroup (Figure 5).

• Download figure
• Open in new tab

Figure 5.

Adjusted (model 1) risk of sustained VT/VF with appropriate ICD therapy associated with sum absolute QRST integral (SAIQRST) in ischemic cardiomyopathy (A; ICM) and nonischemic cardiomyopathy (B; NICM). Restricted cubic spline with 95% confidence interval shows a change in cause-specific Cox hazard ratio (HR) (Y-axis) in response to SAIQRST change (X-axis). The 50^th percentile of the SAIQRST is selected as a reference. Knots of the SAIQRST in ICM patients are at 74-127-186-372 mVms, and in NICM patients are at 80-146-242-462 mVms.

Differences in the risk factors of two competing outcomes

The importance of considering naturally competing outcomes in ICD patients is critically important.²⁶ However, very few studies conducted a formal comparison of the strengths of association of a risk factor with two competing outcomes.²⁷ In our study of primary prevention ICD recipients, unsurprisingly, VT/VF detection duration and VT zone cutoff were distinctly different attributes of the risk of ventricular tachyarrhythmia with appropriate ICD therapies.²⁸ After adjustment for HF characteristics, comorbidities, and management, each additional 15 bpm increase in VT zone threshold was associated with reducing the hazard of sustained ventricular tachyarrhythmias and appropriate ICD therapies by 25%. Long detection duration reduced the hazard of appropriate ICD therapy by 58%. Our results are consistent with the well-recognized fact that appropriate ICD therapies occur more frequently than SCD,²⁹ and ICD programming determines whether patients receive unnecessary therapies for what would have been nonsustained VT with longer ICD detection times.²⁸

We observed that diabetes, BUN, NYHA class, and HF therapy (ACEIs, beta-blockers) were strongly associated with all-cause death without appropriate ICD therapies but did not associate with the primary arrhythmic outcome. These clinical characteristics are often perceived as risk factors of both outcomes. Observed differences in these clinical risk factors’ association with two competing outcomes should be considered in planning future clinical studies and developing SCD risk models.

Age was associated with both outcomes in distinctly different ways. After adjusting for confounders, for every 1 SD (13 years) increase in age, the hazard of death without preceding appropriate ICD therapy increased by 35%, while the hazard of sustained ventricular tachyarrhythmia leading to appropriate ICD therapy decreased by 21%. Age is a ubiquitous risk marker. It is well-recognized that many ICD recipients, despite having an ICD, often die secondary to progressive heart failure or other non-cardiac comorbidities that congregate with advanced age.

Similarities in the risk factors of two competing outcomes

Noticeably, a set of particular clinical risk factors (sex, type of CM, atrial fibrillation, use of AADs, LVEF, eGFR), and ICD type had similar associations with both outcomes, suggesting common or overlapping mechanisms. All of these risk factors are well-known; however, our study conducted a novel formal statistical comparison of the strengths of their associations with competing outcomes in primary prevention ICD recipients. Common mechanisms of two competing outcomes highlight challenges in selecting an ideal candidate for primary prevention ICD.

Global electrical heterogeneity

Our study confirmed that after multivariable adjustment in competing risk analysis, there is an independent association between SVG direction and sustained ventricular tachyarrhythmias leading to appropriate ICD therapies, but not the competing death. The SVG vector points towards the area of the myocardium with the shortest excited state,^{4, 10} and deviations from normal suggest the accumulation of a critical mass of abnormal electrical substrate, which might predispose to ventricular arrhythmias. Our results support the use of GEH, and specifically the orientation of the area SVG vector, as a marker of abnormal underlying electrophysiological substrate responsible for a propensity for sustained ventricular arrhythmias.^{8, 9} Further studies of mechanisms behind SVG directions are warranted.

Notably, we observed that the presence of PVCs on the 12-lead ECG was also associated with the competing risk of appropriate ICD therapy, but not all-cause mortality. Cardiac memory developing in response to PVCs may affect GEH.³⁰ It is well known that PVCs are associated with the development or worsening of HF, and that PVCs can trigger life-threatening ventricular tachyarrhythmias. Further studies of the interaction between an underlying arrhythmogenic substrate (manifesting by GEH) and PVCs (manifesting by cardiac memory development) are needed to uncover the underlying mechanisms, leading to novel therapies.

In contrast to SVG direction, which had significantly different cause-specific hazards associated with two competing outcomes, area SVG magnitude and spatial QRS-T angle had similar associations with both outcomes. The finding of QRS-T angle being associated with both competing outcomes is expected. Numerous previous studies have shown an association between QRS-T angle and broadly defined cardiovascular disease.³¹ Our results underscore the importance and added value of a comprehensive GEH concept, which includes measurement of SVG direction, in addition to the QRS-T angle.

Our results also explain a previously observed discrepancy in SAIQRST. Similar to this study, in another primary prevention ICD study that enrolled NICM patients (PROSE-ICD),^{16, 25} smaller SAIQRST was associated with increased risk of appropriate ICD therapies. In contrast, in the ICM MADIT-II study,²⁴ and the general population,³ larger SAIQRST was associated with ventricular arrhythmias and SCD. The current study demonstrated an opposite association of SAIQRST with sustained VT/VF and appropriate ICD therapies in patients with ICM or NICM.

In our previous studies of healthy individuals¹¹ and the general population,^{7, 22, 23} we primarily observed agreement between findings obtained by using area-based and peak-based GEH metrics. However, in this study of HF patients with reduced LVEF, there were some differences in the results conveyed by peak-based and area-based SVG azimuth and magnitude and QRS-T angle. Theoretically, such differences can be explained by the underlying biophysical and physiological meaning. By measuring peak QRS and peak T vectors, we assess the moment when most of the heart is depolarized (QRS) or repolarized (T).^{32, 33} By measuring QRS and T areas, we aim to assess the complete depolarization (QRS) and repolarization (T) phase.

However, in the diseased/failing hearts, we may not be able to accurately identify an onset and offset of depolarization and repolarization in the entire heart by measuring the onset and offset of QRS and T on the surface ECG due to the heterogeneity and multipolar nature of these processes.³⁴ A single-diploe ECG approximation carries inherent limitation in the modeling of multipolar electrical activity in diseased hearts.³⁴ In the diseased hearts, peak-based measurements can be more meaningful and carry less error than area-based GEH metrics.

Statistical assessment of competing risks

In this study, we preferred to use cause-specific hazards rather than the analysis of subdistribution hazards because Fine and Gray’s method is not valid for causal inference and can produce misleading results.^{35, 36} We used cause-specific hazard that removes an individual from the risk set when any type of event occurs. In contrast, the subdistribution hazard does not remove an individual from the risk set when a competing event occurs. Nonetheless, those who experience one competing event are no longer at risk of the other competing event. As a result, the Fine and Gray method does not isolate distinct causal effects on the competing risks but instead may confound them. The interpretation of Fine and Gray’s method is challenging when outcomes are not rare (as in this study).³⁷ On the other hand, cumulative incidence functions are useful for prediction, which we utilized in our previous study where we developed a GEH-based SCD competing risk.³

Strengths and Limitations

Our study’s strengths include the large sample size and population drawn from a diverse group of medical centers from across the United States. However, because the study included only academic medical centers, the results may not apply to all patients with ICDs. Because of our study’s retrospective nature, ICD programming was not standardized, and differences in ICD detection and treatment parameters could have influenced the results, although we attempted to correct this by including ICD programming in our multivariable models. There was no postmortem ICD interrogation data, and we cannot precisely determine if there were ventricular arrhythmias close to the time of death. Importantly, appropriate ICD therapy is a surrogate forSCD, as all ventricular arrhythmias treated by the ICD would not necessarily have led to cardiac arrest or SCD. As in any observational study, residual confounding cannot be completely eliminated.