VULNERABLE CAROTID PLAQUE: ROBUSTNESS AND CLASSIFICATION CAPABILITIES OF MRI RADIOMIC FEATURES

Patients and image segmentation

Data from fifty-five patients with bilateral carotid stenosis from the Chinese atherosclerotic risk evaluation (CARE-II) multicentre clinical trial were used in this study (ClinicalTrials.gov Identifier: NCT02017756). IRB approvals were obtained for the entire study and for each participating institution, all participants in the study provided written informed consent [7]. The primary objective of CARE-II study was to determine the prevalence and characteristics of high-risk atherosclerotic plaques in the carotid arteries in Chinese patients with recent ischaemic stroke or transient ischaemic attack (TIA). In this study, linear logistic regression models were trained to differentiate culprit and non-culprit carotid artery lesions. The term ‘culprit’ always refers to the causative ipsilateral carotid artery or plaque, as determined following the assessment of ultrasound and MRI examinations. The ‘non-culprit’ side was deemed to have no causative basis for symptoms.

Those recruited underwent MRI of bilateral carotid arteries with the following sequences: three-dimensional (3D) time-of-flight (TOF), T1-weighted (T1) quadruple inversion recovery (QIR), T2-weighted (T2) multi-slice double inversion recovery (MDIR), and Magnetisation Prepared Gradient Recalled Echo (MPRAGE). Detailed MR imaging parameters can be found in the references [7], and in Supplemental Material A. Detailed patient demographics are provided in Table 1 below.

The MRI segmentations were performed manually by experts with reference to 3D time-of-flight (TOF), T1-weighted, T2-weighted and magnetisation-prepared rapid gradient echo (MP-RAGE) images to discern lumen and outer wall boundaries and various atherosclerotic components including IPH, LRNC and calcification [8]. The segmentation and related pre-processing procedures were as follows: (1) lumen and outer wall boundaries were drawn manually on the T1-weighted images to define the region of interest (ROI); (2) all images including the masks for the lumen and outer wall were converted into NIFTI format; (3) radio frequency bias correction [9] and de-noising using a spatially adaptive method [10] were performed; (4) global rigid and local deformable co-registrations were performed with T1-weighted images being the reference modality [11]; and (5) pixels of all images were finally extracted according to the segmentation masks. Except for Step (1), all other steps were performed automatically in MR-VascularView© (Nanjing Jingsan Medical Science and Technology, Ltd., Jiangsu, China).

Radiomic feature extraction and robustness assessment

The PyRadiomics tool (v2.2.0) [12] was used to extract radiomic features based on the MRI pixels from slices which contained atherosclerotic disease defined by local plaque burden ≥ 50%. The feature families obtained by running multi-slice radiomic analysis were: shape features; first-order intensity histogram statistics; Grey Level Co-occurrence Matrix features (GLCM) [13–15]; Grey Level Run Length Matrix (GLRLM) [16]; Grey Level Size Zone Matrix (GLSZM) [17]; Grey Level Dependence Matrix (GLDM) [18]; and Neighbouring Grey Tone Difference Matrix (NGTDM) [19]. The default PyRadiomics setting, a fixed bin width of 25, was used. No resampling was applied.

The overall workflow is illustrated in Figure 1 and includes: calibration; pre-processing procedure to select robust features; development of the predictive model for the identification of culprit lesions; a post-processing procedure to select the most useful features; and the validation methodology. The dependence of radiomic features on the ROI, which defines the area from which features are extracted, was tested by measuring feature robustness using the intra-class correlation co-efficient (ICC) [20–22]. The ROI was manipulated by expanding and eroding the wall boundary outwards and inwards with 1-2 pixels. Only stable radiomic features would then be carried forward to the final analysis and those highly sensitive to minor ROI changes and image quality were omitted.

• Download figure
• Open in new tab

Figure 1.

Schematic which demonstrates the methodology of extracting and evaluating radiomic features, and the cross-validation technique.

Feature regularisation and model cross-validation

Least absolute shrinkage and selection operator (LASSO) and ridge regression were used to reduce over-fitting by introducing regularisation penalty terms to the cost function. In both LASSO and ridge regression, large co-efficients in the linear model are penalised by a characteristic penalty term, though in slightly different ways. LASSO regression involves the L1 norm of the vector of feature co-efficients w in its penalty term, which is the sum of the absolute values of the feature co-efficients. In ridge regression, the penalty term involves the square of the L2 norm, which is the sum of squares of feature co-efficients [23]. The corresponding penalty terms are shown below, where N refers to the number of features involved in the model.

Since the sum of squares of co-efficients less than one is less than their sum, the ridge regression penalty is more accommodating for a more balanced co-efficient vector, while LASSO regression tends to favour multiple co-efficients being set to zero, resulting in a sparser model. While both types of regularisation were performed, with feature selection in mind, the results following LASSO regularisation were typically considered in the final analysis.

To validate the predictive ability of the model and assess its sensitivity to the dataset, a five-fold cross-validation scheme was employed where the data were split into five equally sized folds. Four fifths of the data are used to train a linear logistic regression model, and one fifth of the data are left unseen to evaluate the model. Additionally, to assess the repeatability of this analysis, this scheme was performed on resampled datasets 1,000 times according to the bootstrapping technique of resampling with replacement. Bootstrapping resampling measures the sensitivity to the data set of the cross-validation analysis. This also allowed for a more comprehensive analysis for feature selection.

The combined model and feature selection

The combined model consisted of the radiomic model and conventional risk features composed of morphological and compositional metrics. Morphological metrics included: degree of stenosis; minimum minor axis length (MMAL); minimum lumen area (MLA); inward remodelling index (IRI); and outward remodelling index (ORI), plaque burden, and plaque volume. Compositional features included calcification presence, calcification volume, IPH presence, IPH volume, and LRNC presence, and LRNC volume. The merit of these features alone was considered in the univariate analysis.

To select features, a protocol was developed based on the cross-validation procedure used to identify the most useful features for the combined model. Using LASSO regularisation, the frequency of non-zero feature co-efficients across all five cross-validation folds and all bootstrap iterations was plotted to illustrate the features that were most commonly included across all possible models. With ridge regression, a similar principle was utilised, except that the features were considered selected if their corresponding average co-efficient across all five cross-validation folds was found to be in the top ten among all features considered in the model.

Statistical analysis

Univariate analysis was performed using the χ² test, Fisher’s exact test, or a student t-test when appropriate and those with p-value ≤ 0.1 were included in the multivariate analysis. The ICC was defined using a 2-way mixed analysis of variance (ANOVA) model. The relationship between sensitivity and specificity forms the basis of receiver-operating characteristic (ROC) curves. For all possible cut-offs (the point which separates culprit and non-culprit predictions) of the logistic regression S-curve, a point was added to the ROC curve and the resulting area under the curve (AUC) calculated accordingly. The AUC value is then used to measure the effectiveness of a given predictive model. The significance of the difference between two ROC AUC measurements were assessed using the Wilcoxon statistic for paired samples [24].

Univariate baseline predictors

The baseline model is firstly built upon the consideration of the univariate predictive value of morphological and compositional features. Geometric features except for ORI all had significant univariate predictive value, as did calcification presence and IPH presence. This is presented in Table 2 below. The multivariate predictive value of these features and the significance of the multivariate odds ratio within the model are also presented in Table 2.

Robustness of radiomic features

Robustness of the radiomic features across 12 extraction settings (bin width: 10, 15, 20, 25, 30, 35 and bin count: 8, 16, 32, 64, 128, 256) was calculated for all 107 radiomic features produced by PyRadiomics. The T2-weighted images provided a higher proportion of moderate (0.5 ≤ ICC < 0.9) and high (ICC ≥ 0.9) robust features in comparison with T1-weighted images. Of the features derived from T1-weighted images, 20% had low robustness (ICC < 0.5), compared with almost none which were derived from T2-weighted images. The bin width was a more reliable extraction setting compared to the bin count, which led to a high number of low robustness features. This is presented in Figure 2.

• Download figure
• Open in new tab

Figure 2.

Robustness of T1- and T2-weighted images-derived radiomic features across 6 bin width and bin count settings (Low: intra-class correlation coefficient (ICC)<0.5; Moderate: 0.5≤ICC<0.9; High: ICC≥0.9).

Integration of the radiomics and conventional models

The classification capabilities of models with different features included is tabulated in Table 3 across five cross-validation folds. It can be seen that if only lesion geometry (including stenosis) was considered, the overall accuracy was slightly greater than 60% at the optimal operating point, which is comparable to the case where only lesion components (calcification, etc) were considered.

The combination of lesion geometric and compositional features improved the classification capabilities, albeit without a significant difference when compared with either plaque morphology or composition alone. The radiomics model had better performance than any of the previous models. The combination of radiomic and lesion morphological and compositional features further improved the accuracy to 75.8%. The ROC plots relating to the performance of these models are illustrated in Figure 3.

• Download figure
• Open in new tab

Figure 3.

Receiver operating characteristic curves of the models considering lesion morphological and compositional features, and radiomic features alone and in combination.

The multivariate odds ratios of all features within the combined model are presented in Supplemental Material B, which highlights the relative weight of each feature within the cross-validated predictive models. This contains a conventional model, radiomic model, and a combined model.

In light of prioritising clinical interpretability, a simplified radiomic model was tested, which includes only seven radiomic features: Gray Level Non-Uniformity (GLSZM), Large Dependence High Gray Level Emphasis (GLDM), Small Dependence Low Gray Level Emphasis (GLDM), Difference Variance (GLCM), Informational Measure of Correlation (GLCM), kurtosis (first order), and skewness (first order). The performance of these simplified features is tabulated in Table 4 and the ROC plots may be found within Supplemental Material C, which contrasts the ability of simplified and full models (containing all available features) to distinguish culprit and non-culprit plaques.

Feature selection

Each feature class was analysed individually, and the corresponding ROC curves were created for calculation of the AUC. This is shown in Table 4. The AUC was calculated based on the performance of predictive models on a hold-out set according to the cross-fold validation method. In general, GLDM features derived from T2-weighted images had the best performance, followed by GLSZM and GLCM. NGTDM features had the poorest predictive value, and the poorest robustness profile.

The frequency of selection using LASSO and ridge protocols developed in this study are presented in the Supplemental Materials D and E. First order features skewness and kurtosis appear consistently in feature selection according to the established protocol, despite the first-order feature family in itself not having the strongest predictive performance. The morphological and compositional features including calcification, IPH and plaque burden (PB) were often selected by the linear model under different perturbations of the training data.