Pixel-wise statistical analysis of myocardial injury in STEMI patients with delayed enhancement MRI

Objectives Myocardial injury assessment from delayed enhancement magnetic resonance images is routinely limited to global descriptors such as size and transmurality. Statistical tools from computational anatomy can drastically improve this characterization, and refine the assessment of therapeutic procedures aiming at infarct size reduction. Based on these techniques, we propose a new characterization of myocardial injury up to the pixel resolution. We demonstrate it on the imaging data from the Minimalist Immediate Mechanical Intervention randomized clinical trial (MIMI: NCT01360242), which aimed at comparing immediate and delayed stenting in acute ST-Elevation Myocardial Infarction (STEMI) patients. Methods We analyzed 123 patients from the MIMI trial (62 ± 12 years, 98 male, 65 immediate 58 delayed stenting). Early and late enhancement images were transported onto a common geometry using techniques inspired by statistical atlases, allowing pixel-wise comparisons across population subgroups. A practical visualization of lesion patterns against specific clinical and therapeutic characteristics was also proposed using state-of-the-art dimensionality reduction. Results Infarct patterns were roughly comparable between the two treatments across the whole myocardium. Subtle but significant local differences were observed for the LCX and RCA territories with higher transmurality for delayed stenting at lateral and inferior/inferoseptal locations, respectively (15% and 23% of myocardial locations with a p-value <0.05, mainly in these regions). In contrast, global measurements were comparable for all territories (no statistically significant differences for all-except-one measurements before standardization / for all after standardization), although immediate stenting resulted in more subjects without reperfusion injury. Conclusion Our approach substantially empowers the analysis of lesion patterns with standardized comparisons up to the pixel resolution, and may reveal subtle differences not accessible with global observations. On the MIMI trial data as illustrative case, it confirmed its general conclusions regarding the lack of benefit of delayed stenting, but revealed subgroups differences thanks to the standardized and finer analysis scale.


Introduction
Cardiac magnetic resonance imaging (MRI) has a central role in myocardial infarction experimental and clinical trials, and is frequently recognized as providing the best in vivo surrogate endpoints for infarct size and microvascular obstruction (MVO) (1)(2)(3). The lesion patterns visible on delayed enhancement images vary from ischemia-reperfusion lesions to final scar within the infarct zone on late Gadolinium enhancement (LGE), with or without MVO region on early Gadolinium enhancement (EGE) and LGE and at the acute phase. However, most analyses only assess a tiny part of the imaging data related to these lesions, whose distribution within the myocardium and across slices is complex. They focus on simple global measurements (mainly size, transmurality, and endocardial surface area) (4-6), which drastically limit the quantification of the mechanisms of ischemia-reperfusion and treatment effects. Recent tools (not limited to MRI) offer more detailed assessment and in particular standardized anatomical diagrams (7), but are still insufficiently exploited for the analysis of populations.
Computational tools inspired by statistical atlases can overcome the main bottleneck to a finer analysis of such lesions, namely anatomical differences across subjects and acquisitions (8). First, they allow advanced statistical analysis of myocardial shape differences across a population, as recently demonstrated on acute myocardial infarction patients (9). In addition, as in our study, they allow bringing all the imaging data to a standardized geometry (the "atlas"). This means that much finer statistical comparisons can be performed, up to examining differences in the lesion patterns at each pixel of the myocardium in the EGE and LGE images. This operation is a prerequisite to detect potential inter-subject differences inaccessible and not evaluated by global or regional descriptors of the lesions. In addition, in a logic of personalized medicine, these tools allow understanding and visualizing how a subject is positioned with respect to a given population and/or how she/he responds to a given therapeutic intervention.
In this paper, we re-examined the delayed enhancement images from the Minimalist Immediate Mechanical Intervention (MIMI) trial (10) with such techniques. This multicenter randomized clinical trial aimed at comparing immediate and delayed stenting in acute ST-Elevation Myocardial Infarction (STEMI) patients, two treatment strategies that are still under debate (11,12). By examining commonly-used global descriptors of the lesions, the original MIMI study failed to support delayed stenting against immediate stenting, and even suggested a deleterious effect of delayed stenting on MVO size. However, similar studies reported potential benefits of such intervention regarding infarct size or the incidence of MVO (13,14). Thus, it sounds as a relevant playground to evaluate if a finer analysis, with standardized comparisons up to the pixel resolution, could better characterize infarct and MVO pattern differences across subjects and provide a fresh look at the quantification of the intervention outcome.

Participants
This study retrospectively analyzed the imaging data from the MIMI study (ClinicalTrials ID: NCT01360242), which was a multicenter, prospective, randomized, open-label trial with blinded endpoint evaluation. The study protocol was approved by the local ethics committee (IRB 2010-048), and complied with the Declaration of Helsinki and French laws. All subjects gave written informed consent.
The included patients were adults with symptoms consistent with STEMI ≤12 h, with ST-segment-elevation ≥1 mm in ≥2 contiguous limb leads or >2 mm in ≥2 precordial leads on the electrocardiogram, intended for primary percutaneous coronary intervention. The clinical and angiographic inclusion criteria and the intervention procedure were detailed in Belle et al. (10).
Patients were randomized between two different strategies: immediate stenting or 24-48 h delayed stenting. The primary endpoint was MVO [% of left ventricular (LV) mass] on cardiac magnetic resonance performed 5 days (interquartile range 4-6) after inclusion.

Imaging protocol and labeling
The MRI protocol was performed in multiple centers 3 to 8 days after inclusion on 1.5 T systems. The analysis performed in this paper focused on EGE and LGE images for which LV and lesion segmentation was feasible (Figure 1), obtained according to the detailed MRI protocol provided in the Supplementary Material section of Belle et al. (10).
The myocardial and lesion contours were manually segmented offline by one experienced observer (LP) and controlled by two other experienced observers (MV, PCr) using commercial software (CVI42 v.5.1.0 Circle Cardiovascular Imaging, Calgary, Canada). They consisted of the endocardium, epicardium, infarct (from LGE), and early and late MVO (EGE and LGE, respectively) ( Figure 2). The infarct zone was determined semiautomatically on LGE images using the full-width half-maximum method. MVO corresponded to hypo-enhancement within the infarct zone on the EGE and LGE images, and was segmented manually. The segmented images were exported as Dicom files. Contours were exported in the XML format as.cvi42wsx files. Both images and contours were loaded in Matlab (v.R2016a, MathWorks, Natick, USA) for post-processing and analysis.
The LV-RV junction was manually identified by one experienced observer (ND) using a single landmark on each slice. The myocardial borders around the LV outflow tract were marked by two landmarks to exclude this region on the concerned slices. Finally, the slice locations corresponding to the endocardial apex and the mitral level (the most extreme basal slice with more than 50% myocardium around the blood cavity) were identified. These locations were extrapolated in case the acquisition did not cover such slices.

Standardization onto a common geometry 1
We first parameterized the myocardium by defining radial, circumferential, and long-axis coordinates on each slice, for each individual. They were estimated on images oversampled by a factor 4 to prevent artifacts in case few pixels covered the myocardium, in particular along the radial direction. Their values ranged from 0 to 1, corresponding to endocardium-toepicardium (radial), apex-to-base (long-axis), and the anticlockwise direction starting from the LV-RV junction (circumferential), as illustrated in Figure 3A, top row.
Then, we arbitrarily defined a common reference geometry onto which standardizing the image data, as a semi-ellipsoid with maximal endocardial and epicardial radii of respectively 30 and 50 pixels, represented on 21 slices of 80 × 80 pixels each ( Figure 3A, bottom row).
Finally, images were warped by mapping each patient's coordinates onto the reference coordinates. This was achieved by linear interpolation tailored for data defined on a scattered grid (here, the myocardial coordinates) ( Figure 3B). No extrapolation was performed out of the minimal and maximal slices acquired.

Pixel-wise analysis
Once standardized onto a common geometry, lesions can be compared at each pixel of the myocardium in the EGE and LGE images. We first estimated the representative (average) infarct pattern for each coronary territory and each treatment subgroup, and visualized it at each slice of the reference geometry ( Figure 4A). A Bull's eye representation also summarized this information across slices ( Figure 4B). The Bull's eye consisted of 21 concentric circles corresponding to all the slices of the reference geometry, and 24 segments around the circumference. The following information was displayed: (i) location, estimated as the maximal value of the infarct pattern along the radial direction; (ii) transmurality, estimated as the average value of the infarct pattern along the radial direction; and (iii) statistical differences between the two treatment options, assessed by the pvalue from the Hotelling t-squared test (performed on the pixel map), displayed in a logarithmic color scale.

Link with other characteristics
Despite being standardized onto a common geometry, the infarct patterns still contain a lot of information, which is challenging to visualize and relate to other patient characteristics. We therefore estimated a simplified representation of all the individual subjects in the studied population that can be visualized in two dimensions. To do so, we used the t-SNE (T-distributed Stochastic Neighbor Embedding) algorithm (15), a widespread statistical method for reducing the dimensionality of complex data and visualizing it. Its purpose is similar to the very popular Principal Component Analysis algorithm, except that the latter performs linear dimensionality reduction, which is not fully suited to medical images compared to non-linear dimensionality reduction (16) as performed by t-SNE.
On our data, t-SNE provided a two-dimensional cloud of points where each point corresponds to a sample in the population, while forcing similar patterns to be represented by nearby points, and conversely dissimilar patterns to be represented by distant points. Then, we examined the link between infarct patterns and other variables of interest by coloring each point according to the value Study flow chart. EGE, early Gadolinium enhancement; LGE, late Gadolinium enhancement; STEMI, ST-elevation myocardial infarction. 1 Code and demo data corresponding to this standardization process will be made accessible upon acceptance at: https://github.com/nicolasduchateau/ CMR-imageDataAlignment Duchateau et al. 10.3389/fcvm.2023.1136760  Frontiers in Cardiovascular Medicine of the variable we inspect. We focused on a subset of relevant clinical and therapeutical variables related to patients' risk: coronary artery of concern, ejection fraction, stenting strategy, delay from onset of symptoms to coronary intervention, sex, age, diabetic status, smoking status, and body mass index. This was implemented in Python, using a public open-source version of t-SNE from the scikit-learn library 2 . The retained hyperparameters were two components (for two-dimensional visualization), perplexity of 10 (to define how much neighboring samples are considered for non-linear dimensionality reduction), and Principal Component Analysis initialization (for improved stability).
Before applying t-SNE, infarct patterns were realigned to the left anterior descending (LAD) territory, to better focus on the pattern shape without being confounded by different infarct localizations and analyse all the population at once. This was achieved by rotating all infarct patterns along the circumference such that the average infarct center of a given coronary territory is aligned to the average infarct center for the LAD territory, as previously reported (17).

Global descriptors
We also estimated commonly-used global descriptors of the main lesion characteristics: the proportion of damaged myocardium (for the infarct and the MVO, in %), the average transmurality over the infarcted myocardium (excluding the preserved myocardium, in %), and the endocardial surface area (in %). These characteristics were both estimated before and after standardization to the reference geometry to quantify the consistency of global observations regarding standardization ( Table 1).

Reproducibility analysis
The identification of the LV-RV junction on each slice, and the identification of the apex and mitral valve slice locations were repeated by the same operator and another experienced operator (LP). The reproducibility of these measurements ( Supplementary  Figures S1, S2) and its effect on the statistical analysis of the FIGURE 4 Representative infarct patterns depending on treatment (immediate vs. delayed) and territories (LAD mid or proximal, LCX, and RCA). (A): average pattern for each subgroup, from apex to base (one slice out of two displayed for the sake of clarity). The dashed region corresponds to incomplete myocardium at the basal level. (B): Bull's eye representation that summarizes this information across slices: infarct location and variability, transmurality, and statistical differences between the two treatment options (p-value displayed in a logarithmic color scale). Figure S3) is reported in Supplementary Material.

Statistical analysis
Standard statistical analysis was performed in complement of the pixel-wise analyses exposed in the previous sections, as follows. Continuous variables were expressed as median and interquartile range, with Mann-Whitney U-test used for intergroups comparisons. Categorical variables were expressed in percentage over the number of samples, with Fisher's exact test used for inter-groups comparisons. P-values below 0.05 were considered statistically significant. All the patient and the global lesion characteristics were analyzed using the SPSS software (Version 21.0, IBM Corp, Armonk, NY).

Patient characteristics
The original MIMI study enrolled 160 patients between June 2011 and December 2012, of which 140 patients had interpretable MRI scans (10). In the present study, the datasets of 17 patients were unexploitable with our analysis techniques for technical issues (sub-optimal LV coverage, or image artifacts); the EGE images were not analyzed for these 17 patients and 6 additional patients for which the EGE data were missing. We therefore analyzed the EGE and LGE images for 117 and 123 patients, respectively, corresponding to 65 patients with immediate stenting and 58 patients with delayed stenting (Figure 1). The baseline characteristics of the patients are summarized in Table 2.
The median resolution of the analyzed images was 1.5625 × 1.5625 × 5 mm. The LV ranged over 17 ± 2 slices, with 1 ± 1 extra slices lying out of the acquired stack of slices. Figure 4 displays the representative infarct patterns associated to each territory and treatment. The left side shows these patterns slice by slice from apex to base, while the right side uses Bull's eye views to summarize specific characteristics of the patterns across slices.

Pixel-wise analysis
Infarct patterns (see "location", "transmurality" and "variability" Bull's eye plots for a synthetic view) were roughly comparable between the two treatments for all territories and   across all slices. However, for the left circumflex (LCX) infarcts, higher transmurality was observed with delayed stenting at the lateral locations (average transmurality in these segments: 36% (immediate) vs. 53% (delayed)). A similar trend was observed for the right coronary artery (RCA) infarcts at the inferior/ inferoseptal locations (average transmurality in these segments: 32% (immediate) vs. 38% (delayed)). These subtle observations were confirmed by the "p-value" Bull's eye, which quantified statistical differences between the two treatment groups at each location. The LAD mid and LAD proximal groups respectively had 9% and 13% of myocardial locations with a p-value <0.05, but rather sparsely distributed across the ventricle and not necessarily matching the infarct zone. In contrast, the LCX and RCA groups respectively had 15% and 23% of myocardial locations with a p-value <0.05, including marked statistical differences grouped over a portion of the infarct zone.   Link between infarct patterns and other variables of interest. Each subplot depicts the same cloud of points, which is a simplified two-dimensional visualization (of arbitrary units) obtained with the t-SNE algorithm of how subjects differ regarding their infarct shape: two subjects with similar/ dissimilar infarct shape are placed close/far from each other, as illustrated on subplot (A). Beforehand, infarct patterns of any territory were realigned to match the LAD territory, so that observations focus on the actual pattern shape and not its position around the myocardium. On subplots (B-J), the color code indicates the value of a specific clinical and therapeutic characteristic of the patients. Infarct shapes appear consistent across a given coronary artery territory and differ more between different territories. A slightly similar trend is observed regarding ejection fraction. In contrast, no link is observed for the other variables.  Figure 5 proposes a two-dimensional visualization of the whole population where subjects with similar/dissimilar infarct shapes are mapped close/far from each other, as estimated by the t-SNE dimensionality reduction algorithm. On each subplot, the color code corresponds to the values of a specific clinical or therapeutic variable for each subject, to examine its potential link with the infarct patterns. We remind that to construct this plot, infarct patterns of any territory were realigned to match the LAD territory, so that observations focus on the actual pattern shape and not its position around the myocardium.

Link with other characteristics
We first observe that even when this anatomical alignment was performed, the shapes of infarct patterns showed anatomical consistency across subjects from the same coronary artery. The shapes of LAD infarcts were quite similar, while RCA and LCX infarcts showed more dissimilar shapes that are likely related to broader anatomical and collaterality variations across patients having RCA and LCX culprit arteries. A more subtle grouping of patterns was observed regarding ejection fraction, which might reflect the varying impact of each coronary artery on the cardiac function. In contrast, no specific grouping of patterns was observed for the stenting strategy, confirming the main finding of the MIMI study. The lack of grouping of patterns for all the other selected variables (delay from onset of symptoms to coronary intervention, sex, age, diabetic status, smoking status, body mass index) indicated the lack of bias and the quality of randomization.
Similar observations were obtained by multiple launches of the t-SNE algorithm, which confirmed results against the random state intrinsic to this algorithm. Table 2 examines the five global descriptors of the lesion patterns (area of the infarct, early MVO, and late MVO, transmurality, and endocardial surface area). No statistically significant differences were observed between the two treatments, for all territories, for all-except-one measurements before standardization / for all after standardization to the common reference geometry.

Reproducibility analysis
Supplementary Figure S1 summarizes the intra-and interoperator variability in the identification of the LV-RV junction on each slice, for the EGE and LGE images. Angles were defined by the lines joining the measurement of each operator and the center of the LV cavity, on each slice. Intra-operator differences were low, with higher differences near the apex but variability in the range [−20.0°,8.7°]. Inter-operator differences were higher but moderate, with higher differences near the apex but variability in the range [−18.9°,23.4°]. Observations were comparable for the EGE and LGE images.
Supplementary Figure S2 complements this analysis by summarizing the intra-and inter-operator variability in the identification of the apical and basal slices, for the EGE and LGE images. Differences were low, mostly ±1 slice both for the EGE and LGE images. They were higher for the intra-operator experiment, in particular at the basal level (more cases with a difference of ±2 slices, a few up to ±3 slices), in part due to the difficulty of extrapolating this location when the LV was not fully covered by the slices.
Finally, Supplementary Figure S3 provides a display similar to the Bull's eyes from Figure 4 for the location, transmurality, and variability of the infarct patterns, and for the p-value quantifying statistical differences between the two treatment groups, for identification of the LV-RV junction and basal/apical slices repeated by the same operator or by a second one. Patterns are very similar to the ones in Figure 4, and lead to equal conclusions: comparable patterns between the two treatments for all territories and across all slices, but higher transmurality with delayed stenting for the LCX and RCA infarcts, confirmed by statistically significant differences near the infarct zone, contrary to the LAD mid and LAD proximal infarcts.

Discussion
We proposed a fresh look at the characterization of myocardial lesions (infarct and early/late MVO) to go beyond global assessment with quantitative and standardized comparisons up to the pixel resolution. On the delayed enhancement imaging data from the MIMI study, which aimed at assessing acute myocardial infarction patients while comparing different reperfusion strategies, our results not only confirmed at a much finer resolution and in a more standardized manner that immediate and delayed stenting led to comparable infarct and reperfusion lesions, but revealed a deleterious effect of delayed stenting at specific locations. They also indicated that the infarct shape was not related to other clinical and therapeutic variables of interest.
Comparing delayed vs. immediate stenting highly depends on the chosen MRI endpoints (2). In the original MIMI study (10), the primary endpoint was MVO (% of total LV mass observed by MRI 5 days after inclusion). Primary endpoints were slightly different in similar randomized controlled trials of moderate sample sizes. The INNOVATION study used the infarct size (% of total LV mass observed by MRI 30 days after inclusion) (14), while the DEFER-STEMI study considered the incidence of no-/ slow-reflow (Thrombolysis In Myocardial Infarction ≤2, 2 days after inclusion) (13). They concluded that delayed stenting may better preserve the myocardium, not necessarily for all territories or supported by statistically significant differences. The much larger DEFER-DAMANI randomized controlled trial did not find evidence of improvement with delayed stenting using a composite score of clinical outcomes (18). Its findings were confirmed by a sub-study, which considered the final infarct size as primary endpoint (19). Current guidelines do not therefore Duchateau et al. 10.3389/fcvm.2023.1136760 Frontiers in Cardiovascular Medicine recommend delayed stenting (12). Nonetheless, differences in the conclusions brought by these studies should be tempered, mainly due to different primary endpoints, different delays in the stenting procedure, and different spectrum of included patients (11). Delayed stenting might still benefit a subset of patients, to be confirmed by other randomized controlled trials and reanalysis of the existing studies as targeted in the PRIMACY trial (20).
In the literature, the lack of clear differences between treatment options may also come from the simplicity of global descriptors (occurrence, size, transmurality, or endocardial surface area as a surrogate of the area at risk, mainly) to assess lesions of complex shapes. Many studies investigating the cardiac function, not necessarily with MRI, already underlined the value of regional or local descriptors of the observed diseases (7,16). Our approach also goes into this direction and provides a more complete and standardized assessment of the lesion patterns up to each pixel of the myocardium. In our study, this local quantitative assessment not only confirmed the lack of benefit from delayed stenting, but also revealed subtle differences insufficiently rendered by global measurements and suggesting a deleterious effect visible at some locations (in particular, higher transmurality with delayed stenting for the LCX and RCA).
Statistical atlases have been widely used to quantify shape differences across of a population (8), as recently demonstrated on acute myocardial infarction patients (9). In our work, we used these techniques beyond shape assessment and analyzed the imaging data on a common geometry. This standardization was inspired by previous works on the statistical analysis of myocardial infarct (21,22), except that these considered 3D mesh data (while we directly operated on the image data), and did not incorporate MVO. The use of such spatial alignment goes in the sense of recent recommendations, to better standardize the imaging data (2).

Study limitations
Seventeen patients were not analyzed due to technical issues related to the image quality. Besides, our population size was moderate, and the statistical analysis of lesion patterns may be affected by this sample size, in particular in case of small or unbalanced subgroups as observed for some territories. We therefore supported as much as possible our observations with complementary measurements (the variability and statistical differences computed locally, the global lesion characteristics, and intra/inter-operator reproducibility analysis).
Spatial alignment required the semi-automatic segmentation of the myocardium and the manual identification of few landmarks on each slice, which is time-consuming and needs careful quality control in the current implementation. This could be solved by considering recent advances in the automatic segmentation of delayed enhancement MRI, when its accuracy reaches acceptable ranges on all slices on routine imaging data (23,24).
Non-linear statistical analysis techniques may overcome some limitations of the techniques we used here to estimate representative lesion patterns (21,25), but would require larger populations and may still provide limited rendering of sharp transitions in the analyzed data, such as the interface between the myocardium and the infarct (17).

Conclusion
Our methodology allows assessing the acute myocardial infarction lesions much beyond commonly used global descriptors while delivering a synthetic picture of a population, and could be used for standardized reporting in many cardiac MRI clinical research studies. This approach provides quantitative finer insights into how lesions differ between subgroups of subjects, and enables to locate a given subject within the population/group response in a logic of personalized medicine. On the MIMI study, it confirmed at a much finer scale the comparable myocardial damage between immediate and delayed stenting, and even suggested a deleterious effect of delayed stenting visible at some locations. This strategy appears highly promising to analyze in a much more integrated manner the multi-parametric and longitudinal data from ischemia-reperfusion studies, and better state on potential therapies to reduce myocardial loss after myocardial infarction.

Data availability statement
Data generated or analyzed during the study are not publicly available. Requests to access these datasets should be directed to nicolas.duchateau@creatis.insa-lyon.fr.

Ethics statement
The studies involving human participants were reviewed and approved by the local ethics committee (IRB 2010-048). The patients/participants provided their written informed consent to participate in this study.
French ANR (MIC-MAC project, ANR−19-CE45-0005), the Fédération Francaise de Cardiologie (MI-MIX project, Allocation René Foudon), and the Institut Universitaire de France. These funding sources had no involvement in study design; in the collection, analysis and interpretation of data; in the writing of the report; nor in the decision to submit the article for publication.