Cardiotoxicity detection tool for breast cancer chemotherapy: a retrospective study

Code development

In this study, we developed a novel package of codes, the “RNA Toolbox” to derive potentially useful parameters of cardiotoxicity using semi-automated techniques to enhance reliability. MATLAB was employed to develop codes performing core functions such as filtering, correction, and segmentation (using phase and amplitude images), to extract features from the RNA images.

Dataset

We retrospectively analyzed 110 patients with breast cancer who were treated with various cardiotoxic agents. Patients younger than 75 years with good RNA image quality, no missing clinical information, and who were not terminally ill from breast cancer were included in the study. Patients with cardiac illnesses (n = 1), with breast cancer beyond stage 3 (n = 3) or underwent radiotherapy (n = 3) were excluded from the study. Seventy-seven patients underwent cardiotoxic 5 Fluorouracil (5FU), Epirubicin, and Cyclophosphamide (FEC) chemotherapy, while the remaining patients received a combination of trastuzumab biological therapy and FEC therapy at Aberdeen Royal Infirmary, NHS Grampian, Aberdeen, UK. All patients underwent a baseline scan, and one or more RNA scans were performed after chemotherapy. All RNA images were assessed using a validated clinical package called eSoft^® (Siemens Healthcare GmbH, Erlangen, Germany), along with the in-house MATLAB^® RNA Toolbox (The Math Works, Inc. Natick, NY, USA; Version 2020a). This study was approved by the Institutional Ethics Review Board of the University of Aberdeen and Aberdeen Royal Infirmary, NHS Grampian. The requirement for informed consent was waived due to the retrospective nature of this study.

Optimization process

The RNA Toolbox requires optimization regarding background subtraction, image correction, filtering, and image segmentation.

Background optimization

Background subtraction corrects noisy RNA images, but always requires optimization. To analyze the effect of background spatial location, we wrote a MATLAB code that automatically drew six background regions of interest (ROI) in six different spatial regions outside the heart area (Fig. S1). Thirty-two images were analyzed to define the best region from which to draw the background ROI and further optimize the RNA Toolbox. After multiple trials, a background width of four pixels was observed to produce the mean gray values most representative of the background around the heart. Therefore, the background width was standardized to four pixels in all scenarios. We compared the background mean gray values with each other and with those produced from clinical software (eSoft^®) using the Pearson correlation. The background in eSoft^® was automatically drawn close to the distal boundaries of the LV (location 1 in Fig. S1). A repeated-measures analysis of variance (ANOVA) was performed to compare the mean gray values of the corresponding backgrounds to account for the random effect of the individual.

Filter optimization

This step aimed to evaluate the effects of different filters on the three parameters, LVEF, E, and ApEn, to optimize image processing and filtering using the RNA Toolbox. These parameters were selected to evaluate the performance of the segmentation techniques for each filter scenario. This assesses the filtering within the segmentation context. Four filters (median, mean, Wiener, and median-modified Wiener filter (MMWF)) were compared under two different simulations: 3 × 3 and 5 × 5 kernels. For this purpose, repeated-measures ANOVA and intraclass correlations were used to compare the different filters. Finally, a repeated-measures ANOVA was performed to evaluate the effects of changing filter type and corresponding kernel size on E, ApEn, and LVEF. The eSoft^® filter was set to a Butterworth filter with a cut-off of 0.4 and an order of 5.

Optimization of segmentation technique

We developed a novel technique to accurately segment the LV in an RNA scan consisting of 16 sequential frames representing the mean heart volume over time. LVEF was measured using custom MATLAB^® code utilizing the novel segmentation algorithms (ACRG). The ACRG used phase and amplitude images because these provide better contrast than spatial domain RNA images. The steps utilized in the ACRG for LV segmentation are outlined below.

An initial seed point was manually selected as the starting point for calculating the initial contour using the RG algorithm. This produces 16 initial contours over 16 frames, which typically extend outside the LV. Using these contours and depending on the grey values within them in each frame, the AC algorithm begins to work on the phase image.

In the phase image, using the AC algorithm, the initial contour evolves over time (in a shrinking pattern) until it finally reaches the LV boundaries in each frame. A set of new contours (i.e., 16 masks) was generated for the next step.

The new contours were applied to the amplitude image to define the maximum boundaries for the RG algorithm. Subsequently, the algorithm began to evolve using the same original seed point from the first step.

Following this, even more precise contours were generated and applied to the actual image to calculate volumes and corresponding activities. Step 2c in Fig. 1 illustrates this occurrence.

LVEFs derived using this novel segmentation algorithm were compared with readings from clinically validated eSoft^® software using appropriate methods to address consistency and reproducibility. Agreement between the eSoft^® and MATLAB^® LVEFs was evaluated using three numerically based processes: intraclass correlation coefficient (ICC), Bland & Altman diagram, and Lin’s concordance correlation coefficient (LCCC).

Predictive models

Continuous variables were expressed as mean ± standard deviation for normally distributed data and as median with interquartile range for non-normally distributed data. Univariate analysis was performed using linear regression to assess the association of outcome (LVEF) with the following predictors: E, ApEn, bounded-ApEn, S, LHR, circularity (diastolic and systolic), elongation (elongation of diastole and systole), short/long-axis FS (diastolic and systolic FS), mean beat duration, rejected/accepted beats, heart rate, or phase angle. Factors associated with LVEF (p < 0.2) in unadjusted univariate linear regression were included in a stepwise multivariate model to identify independent factors associated with LVEF, and the variables were assessed for their usefulness in predicting LVEF. The coefficient of determination (R²) was calculated to determine the proportion of variation in LVEF explained by a set of independent variables Depending on the data, an appropriate post-hoc analysis was conducted to determine the significance of each levels of a variable and the p-values from a test were adjusted by employing Bonferroni correction to account for multiple comparisons. Estimates of beta coefficients and 95% confidence intervals (CIs) were evaluated to determine the extent of association for each independent variable. Linearity and homoscedasticity were examined using scatter and Q-Q of standardized residual plots, as well as Kolmogorov–Smirnov and Shapiro–Wilk tests. The multicollinearity between predictor variables was evaluated using the variance inflation factor (VIF), with values greater than 10 indicating its presence. A multiple linear regression model was used to evaluate the relationship between post-chemotherapy LVEF and all other predictors, derived from the pre-chemotherapy RNA scans. Additionally, the novel RNA Toolbox was used to analyze 103 RNA scans to evaluate the receiver operator characteristic (ROC) curves. See Fig. 1 for further insights into the methods used in this study. All analyses were conducted using SPSS for Windows (version 27; Chicago, IL, USA).

Patient characteristics

The median age of the 103 patients included in this study was 37 (20–81) years. All patients had either stage 2 (n = 79; 39%) or stage 3 (n = 24; 61%) breast cancer with no pre-therapy history of cardiac abnormalities. The corresponding chemotherapy dosages are summarized in Table S2.

Image background optimization

Image backgrounds (Bkg) 1–6 were automatically drawn around the heart along the axis, as shown in Fig. S1. The mean gray values between at least three groups (F (3.5, 415) =10.3, P < 0.001) were statistically significant. Pairwise comparisons revealed that the difference in the mean gray value counts between Bkg1 and Bkg2 (P < 0.05), Bkg1 vs. Bkg6 (P < 0.05), and Bkg2 vs. Bkg3 (P < 0.001) were statistically significant (Table 2). Table S3 shows mean gray values of the corresponding backgrounds with standard deviations. When comparing experimental backgrounds (i.e., Bkgs 1–6) with the eSoft^® Bkg, the lowest P-value was found between eSoft^® Bkg and Bkg2. No significant differences were observed between eSoft^® Bkg and Bkg1, Bkg3, Bkg4, Bkg5, or Bkg6 (P = 0.98). Overall, the results suggest that a change in background spatial locations substantially affects mean gray value. Pearson’s correlation was used to analyze the associations between different backgrounds. The inter-background correlation matrix depicts various correlation values between the ROIs, indicating that changing the background ROI location may affect mean gray values. The strongest correlation was observed between the eSoft^® background ROI and Bkg1 (Table S4).

Filter optimization: a comparison of median, mean, Wiener, and MMWF filters

Our null hypothesis was that the distributions of the median, MMWF, Wiener, mean, median3, MMWF3, Wiener3, and mean3 would be the same. Filter names without numbers (e.g., median filter) represent filters with a $5\times 5$kernel size, whereas filter names with “3” as a suffix (e.g., median3 filter) represent filters with a $3\times 3$ kernel size. Friedman’s test of differences among repeated measures was conducted, and a chi-squared (X²) value of 28 were obtained. A statistically significant difference was observed in the LVEF distributions with different filters (X²(7) = 28; P = 0.001).

Post-hoc analysis using Wilcoxon signed-rank tests was conducted, with the Bonferroni correction applied at all significance levels. The descriptive statistics for all filters and the corresponding kernel sizes are summarized in Table 3. Significant differences were observed between the mean3 and median filters (Z = 2.11, P = 0.002) and between the Wiener3 and median filters (Z = −2.08, P = 0.002). Despite an overall change in LVEF, no significant difference in the LVEF values was observed. Moreover, all significant differences appeared to arise between kernel sizes rather than filter types. Thus, post-hoc tests were performed individually for each kernel size. The results showed no significant differences between the individual 5 × 5 and 3 × 3 kernel size filters (P < 0.05), suggesting no evidence of a difference in the corresponding LVEF distributions of the filters when kernel parameters were fixed. It may suggest an association of LVEF with kernel size but not with filter type. Spearman’s correlation test was performed between $3\times 3$ and $5\times 5$ filter vs. eSoft^® LVEFs. The highest correlation between eSoft^® and all experimental LVEFs was recorded for the 5 × 5 median filter (r = 0.82).

A repeated-measurement ANOVA was performed to evaluate the effect of collectively changing the filter type and corresponding kernel size on entropy (i.e., the interaction effect). As repeated-measures ANOVA assumes independence, normality, and sphericity of the data, these assumptions were tested. The Kolmogorov–Smirnov and Shapiro–Wilk normality tests showed the normal distribution of all variables. Mauchly’s test results showed that our E data did not meet the sphericity assumption (P < 0.001) (Table 4). The Greenhouse–Geisser results identified a statistically significant effect of filter (F(1.2, 126) = 288.446, P < 0.001), kernel size (F(1,42) = 81.147, P < 0.001), and filter × kernel interaction (F(1.953,126) = 126.540, P < 0.001) (Table 5). To follow up on these significant main effects and consider the significant interaction (i.e., filter × kernel size), a paired sample t-test was performed for each filter with the corresponding kernel size. The paired sample statistics are presented in Table 6. A greater difference in the mean E values was observed for the median and mean filters (when changing the kernel size) than for the other filters. Thus, changing the kernel size in the median, mean, and Wiener filters significantly affected the E values (P < 0.001) (Table 7). Both kernel sizes depicted alternating patterns in all filters, revealing a significant difference in the mean entropy values (E) between the filters, except the MMWF (Fig. 2). Thus, the MMWF ( $5\times 5$ kernel) was selected for the remainder of this study.

Optimizing segmentation: Bland-Altman test (eSoft^® vs. ACRG)

The graphs in Fig. 3 show the LVEF obtained from 43 patients using two different methods: eSoft^® and the RNA Toolbox using ACRG under four different filters. Under all filter settings, a funnel effect was observed, in which differences were reduced with an increase in LVEF values. Additionally, the mean difference deviated significantly from zero, indicating that the RNA Toolbox ACRG provided systematically higher values than eSoft^® for LVEFs (see Table 8). The funnel effect does not indicate disagreement, but rather a distribution effect; in addition, this effect may have occurred due to the limited sample size. Thus, LCCC and ICC tests were performed to check for correlations.

Optimizing segmentation: LCCC (eSoft^® vs. ACRG)

The estimated slope of the regression line (eSoft^® vs. ACRG) passes through the center of data at 0.98 (under all filter settings), suggesting high correlation of eSoft^® values and RNA Toolbox ACRG outputs. In addition, an average value of R² = 0.90 suggested that 90% of the variation in the eSoft^® variable could be explained by its linear relationship with the RNA Toolbox ACRG. The average value of the concordance correlation coefficient (0.90) was also close to 1 (Table 9), showing good agreement between the two methods. Moreover, the heatmap in Fig. 4 shows hot spots dispersed and located along the line of best fit, indicating that this method may have good agreement with eSoft^®.

Optimizing segmentation: ICC (eSoft^® vs. ACRG method)

The reliability of a single rating was estimated (number of scans = 43). The estimated mean value of Pearson’s correlation coefficient of all LVEF readings under different filter settings was 0.9 with a 95% CI [0.97–0.98], which was almost identical to the corresponding LCCC results. As the ICC formula considers systematic effects, an ICC value of 0.95, which is very close to Pearson’s correlation coefficient, was obtained. The relatively high ICC may represent a good agreement between the two methods. Figure 5 shows multiple variables graphs of eSoft^®-derived LVEFs vs. RNA Toolbox ACRG-derived LVEFs under the four filter settings.

Regression model

The results of the general regression model (n = 103) showed that 41% of the variance in LVEF could be accounted for by five predictors, collectively, (F (5, 97) = 15, P < 0.001). Looking at the unique individual contribution of the predictors, the results showed that LVEF (derived from phase image) $(\beta =0.09,P=0.001)$, LV ejection rate $(\beta =4,P=0.03)$, ApEn $(\beta =11.6,P=0.001)$, ApEn (D&S) $(\beta =39,P=0.002)$and LV systole size $(\beta =0.03,P=0.01)$ were significantly associated with LVEF. These predictors reported acceptable VIF (<1.4), suggesting no serious multicollinearity concerns. The residual plot of this model shows a relatively symmetrical distribution with a tendency to cluster toward the middle of the plot in a random pattern, indicating a good fit for the linear model (Fig. 6). This suggests that cardiotoxicity resulting from anthracycline chemotherapy in patients with breast cancer may affect LV volume, size, and LV dimensions.

The analysis of results indicates that count-based and size-based predictors made a significant predictive contribution in the model, with LVEF being notably influential in the model. Time-based predictors, such as rejected beats and mean beat duration have less contribution to the model and were characterized as non-image-based parameters. These criteria were determined based on the inclusion of normal/abnormal values and their respective sources.

Risk stratification for chemotherapy-related cardiotoxicity

All predictors were implemented in an ROC curve, which is a test for evaluating the diagnostic performance of predictors and indicates how much a model is capable of distinguishing between classes We re-evaluated the model’s performance in predicting abnormal LVEF (<54%) (adopted from literature (Hesse et al., 2008)) after chemotherapy and found it to be moderately successful. The ROC curve results showed that the model was able to distinguish between the two groups, that is, normal (>54%) and abnormal (<54%) LVEF) at a moderate level (see Fig. 7 for the corresponding area under the curve (AUC) value). Furthermore, the phase LVEF, LA FS, and LV ejection rate were shown to significantly contribute to the prediction (P < 0.05).

Limitations

This study was principally limited by its retrospective nature, which makes it prone to misclassification bias. Additional selection bias might have been present as our controls were primarily recruited through convenience sampling, and hence may not represent the general population. The limited sample size and limited accessibility for clinical data were also limitations; therefore, prospective studies are required. This study considered only two scans, which may have exaggerated any bias in the data. The average duration of chemotherapy or biological treatment was 12 months, with a substantial variation in the treatment dosage. However, cardiotoxicity may have occurred after 12 months, and significant data may have been missed. Moreover, the effect of comorbid factors on LVEF reduction following chemotherapy has not been evaluated thoroughly owing to difficulties in accessing a full patient history. Although our segmentation method shows good agreement with the clinical eSoft^® software, a larger sample number is recommended for further validation.