Subjects

We retrospectively identified 26 subjects (12 females, 14 males), referred for clinical CMR imaging between November 2013 and July 2020 for known or suspected HCM, for whom a comprehensive MRI study including both T1 and T2 mapping sequences was performed. The HCM diagnoses followed the latest European Society of Cardiology guidelines and were based on the presence of LV wall thickness ≥ 15 mm in one or more myocardial segments, not explained solely by loading conditions². Clinical characteristics of the patient group are detailed in Table 1. The study was approved by the local ethics committee of the Azienda USL Toscana Nord Ovest (Pisa, Italy). Written informed consent was obtained from all subjects, and all experiments were performed in accordance with relevant guidelines and regulations.

CMR imaging

All CMR imaging examinations were performed using a 1.5 T MRI scanner system (MAGNETOM Avanto, Siemens Healthcare, Erlangen, Germany) equipped with 45 mT/m gradients strength and a 12-channel phased array coil.

Cine images were performed in the 2- and 4-chamber view planes (3 slices each) and in short-axis view (8–14 slices encompassing the entire LV), using a TrueFISP sequence (TR = 2.5 ms, TE = 1.2 ms, slice thickness = 8 mm).

Both T1 and T2 maps were obtained in the short axis view (a single slice located where myocardial thickness, evaluated with cine images, was maximum). T1 mapping was performed utilizing a modified look-locker inversion recovery (MOLLI) pulse sequence with a 3–3–5 acquisition scheme⁶¹. Pulse sequence parameters were as follows: TE/TR = 1.14/2.5 ms, flip angle = 35°, matrix size = 126 × 192, in-plane resolution ranged from 1.77 mm × 1.77 mm to 2.34 mm × 2.34 mm, typical field of view = 380 mm × 275 mm, slice thickness = 8 mm. T2 maps were obtained utilizing a T2-prepared TrueFISP sequence⁶² with the following parameters: T2 preparation time = 0/24/55 ms, TR = 4 × R–R, flip angle = 70°, matrix size = 126 × 192, in-plane resolution ranged from 1.77 × 1.77 mm to 2.34 × 2.34 mm, typical field of view = 380 mm × 275 mm, slice thickness = 6 mm.

Then, 10–15 min after intravenous administration of 0.2 mmol/kg of gadolinium DTPA (Magnevist, Schering), gadoteric acid (Dotarem, Guerbet) or gadoteridol (Prohance, Bracco), LGE images were acquired in 2- and 4-chamber view of the LV (3 slices each), as well as in the short-axis view (8–14 slices encompassing the entire LV) using a 2D phase-sensitive inversion recovery (PSIR) sequence (TR = 700 ms, TE = 1.09 ms, slice thickness = 8 mm, inversion time = 200–300 ms, typical in-plane resolution = 2.4 mm × 3.2 mm).

Preprocessing of T1 and T2 maps

For each subject, a region of interest (ROI) covering the entire myocardium was manually segmented by a single expert radiologist in cardiac MRI, with 15 years of experience, using 3D Slicer (Version 4.11.2)^63,64. ROIs were delineated on T1 and T2 maps separately, avoiding voxels with potential partial volume effect. Supplementary Figure S1 shows an example myocardium ROI for a representative HCM patient.

Given that T1 and T2 maps allowed acquiring only a single slice, resampling voxel size was performed through 2D interpolation by using the B-spline interpolation algorithm (with the origins of interpolation and original image grids aligned together⁵⁹). The original in-plane spatial resolution ranged from 1.77 mm × 1.77 mm to 2.34 mm × 2.34 mm across patients for both T1 and T2 maps. Therefore, T1 and T2 maps were resampled to an in-plane isotropic spatial resolution of 1.8 mm, 1.9 mm, 2.0 mm, 2.1 mm, 2.2 mm, 2.3 mm, and 2.4 mm. Resampling voxel size was performed for different bin widths (see below).

T1 and T2 maps discretization was carried out following a fixed bin width approach, as recommended by the Image Biomarker Standardisation Initiative (IBSI) when dealing with quantitative data⁵⁹. In particular, for each resampling voxel size (i.e., 1.8, 1.9, 2, 2.1, 2.2, 2.3, and 2.4 mm in-plane isotropic spatial resolution), bin width values of 3.60 ms, 3.95 ms, 4.30 ms, 4.65 ms, 5.00 ms, 5.35 ms, 5.70 ms, 6.05 ms, and 6.40 ms were employed for T1 maps, while bin width values of 0.49 ms, 0.50 ms, 0.51 ms, 0.52 ms, 0.53 ms, 0.54 ms, 0.55 ms, 0.56 ms, and 0.57 ms were employed for T2 maps. Indeed, given that the median (across HCM patients) range of T1 and T2 values (i.e., the difference between the maximum and minimum T1 and T2 values of voxels within the ROI) was 265 ms and 26 ms, respectively, bin width values were chosen in such a way that the number of quantization levels of T1 and T2 maps was within the range of 30–130, for each HCM patient. This approach, adopted also in previous technical studies^40,41,52,54, has the potential to limit variability in radiomic features estimation^43,65,66.

Different spatial filters were applied on T1 and T2 maps: 2D wavelets (Daubechies 3) yielding four filtered maps (i.e., wavelet-LH, -HL, -HH, -LL, where L/H refers to the combination of low-/high-pass filters applied in the horizontal and vertical direction), gradient magnitude of the map (i.e., gradient filter), the square of the map values (i.e., square filter), and the square root of the absolute map value (i.e., square-root filter). Specifically, in order to restrict the large number of possible combinations, image filtering was performed at fixed resampling in-plane isotropic spatial resolution of 2.1 mm, and at fixed bin width of 6 ms and 0.56 ms for T1 and T2 maps, respectively, yielding median (across subjects) numbers of quantization levels between 30 and 130^43,65,66.

All preprocessing steps (i.e., T1/T2 maps resampling, discretization, and spatial filtering) and subsequent radiomic features estimation were carried out by using the open source PyRadiomics library⁶⁷ (Version 3.0.1) with Python (Version 3.7.3), running on a MacBook Air (macOS Version 10.14) with a 1.8 GHz Intel Core i5 CPU.

Radiomic features estimation

Given that the used acquisition sequences allowed obtaining T1 and T2 maps on a single slice, the 2D versions of radiomic features were considered. For each ROI and preprocessing combination, in terms of resampling voxel size and bin width, a total of 98 features were obtained: 9 2D-shape features, 16 first order features (14 intensity-based statistical features and 2 intensity histogram features, namely Entropy and Uniformity), and 73 second order features (i.e., textural features) from gray level co-occurrence matrix (GLCM, 22 features), gray level run length matrix (GLRLM, 16 features), gray level size zone matrix (GLSZM, 16 features), gray level dependence matrix (GLDM, 14 features, with coarseness parameter α = 0), and neighborhood gray tone difference matrix (NGTDM, 5 features). Second order features estimation was performed according to the Chebyshev norm with a distance of 1 pixel. GLCM and GLRLM features were computed from each 2D directional matrix (i.e., at 0°, 45°, 90°, and 135°) and averaged over 2D directions.

For each ROI and spatial filter applied on T1 and T2 maps, a total of 89 features were estimated, i.e. all the above features but the shape features. Indeed, given that shape features are usually estimated regardless of the applied image filter, they were not included in our analysis.

All radiomic features were computed in accordance with the definitions provided by the IBSI, with shape features computed in 2D instead of the proposed 3D version. It is worth noting that the first order feature of Kurtosis calculated by PyRadiomics was in accordance with the IBSI except for an offset value (i.e., 3).

Statistical analysis

In this study, three different effects on radiomic features estimation were assessed for T1 and T2 maps singularly:

1. A.

   For each bin width, the effect of using different resampling voxel sizes.

2. B.

   For each voxel size, the effect of using different bin widths.

3. C.

   At fixed resampling voxel size and bin width, the effect of using different spatial filters.

For each effect of interest, any variability in radiomic features estimate was assessed through ICC analysis^68,69. In particular, the two-way mixed effects model, with single rater and absolute agreement, was selected. Accordingly, the ICC coefficient was calculated as:

where MS_R = mean square for rows (i.e., between subjects), MS_E = mean square for error, MS_C = mean square for columns (i.e., between measurements), n = number of subjects and k = number of raters, with ICC matrices realized considering each resampling voxel size, quantization bin width or image filter as a rater. ICC values range between 0 (maximum variability) and 1 (minimum variability) and express the variability of radiomic features estimate associated with the effect of interest with respect to the variance between subjects (i.e., a relative variability). Radiomic features were stratified based on their degree of relative variability^37,52: high (ICC ≤ 0.5), considerable (0.5 < ICC ≤ 0.75), moderate (0.75 < ICC ≤ 0.9), and low (0.9 < ICC ≤ 1) relative variability.

In order to characterize in greater detail and absolutely the variability in radiomic feature estimate, an additional analysis of the coefficient of variation (CV) was performed. In particular, CVs were computed as the percentage ratio between standard deviation and absolute mean values of features estimate obtained by varying one element of interest (resampling voxel size, bin width or spatial filter) while fixing the others. Accordingly, this yielded 26 × 9 CVs for effect A (i.e., 26 subjects and 9 fixed bin widths), 26 × 7 CVs for effect B (i.e., 26 subjects and 7 fixed resampling voxel size), and 26 CVs for effect C (i.e., 26 subjects with fixed resampling voxel size and bin width).

Any linear correlation between radiomic features estimate and resampling voxel size or bin width was assessed by a repeated measures correlation analysis, namely rmcorr⁷⁰. This technique accounts for non-independence among observations (i.e., repeated measures of the same radiomic features of the same subject by varying one preprocessing element) using analysis of covariance (ANCOVA) to statistically adjust for the variability between subjects. A p-value < 0.05, adjusted for multiple comparisons using Bonferroni correction, was considered significant.

Statistical analysis was carried out by using R (Version 3.6.2) software package in the RStudio (Version 1.2.5033) environment⁷¹.