Capabilities of 3D Mechanics and Radiomics Analysis in Predicting High-Risk Carotid Plaques

1. Introduction

Carotid disease is a major subset of cardiovascular diseases (CVD), which remains one of the leading causes of mortality worldwide (World Health Organization, 2021). The current gold standard for risk stratification, luminal stenosis, is typically assessed through linear measurements at the narrowest lumen and referenced to either the normal (NASCET) or estimated (ECST) vessel caliber (Mallarini and Saba, 2010). Patients with severe stenosis (≥70%) often undergo carotid endarterectomy (CEA) or carotid artery stenting (CAS). However, clinical events are more frequently associated with mild-to-moderate stenosis (<70%), where the benefits of surgical intervention are less clear (Naylor, 2008). This highlights an urgent need for advanced diagnostic methods capable of evaluating plaque vulnerability beyond the limitations of stenosis-based assessments.

Advancements in functional magnetic resonance imaging (MRI) have facilitated carotid plaque analyses, particularly in identifying inflammation-related biomarkers (Howarth et al., 2009; Sadat et al., 2011; Tang et al., 2009). Imaging features such as ultrasound-measured plaque area (Nicolaides, 2010) and MRI-detected intra-luminal thrombus (ILT) presence (Gupta et al., 2013) have been consistently recommended for identifying patients at heightened risk of stroke during best medical treatment (BMT) (Aboyans et al., 2017; Naylor et al., 2023). These imaging biomarkers may serve to refine surgical decision-making alongside patient-specific factors such as life expectancy and medication history in managing symptomatic and asymptomatic carotid artery disease.

This study introduces an integrative approach to carotid plaque vulnerability assessment, combining radiomics-based image texture analysis with mechanical modeling of the vessel wall. Radiomics is a data-driven technique that analyses imaging-derived features, such as pixel intensity distribution and heterogeneity, to characterize lesion phenotypes (Gillies, Hricak, and Kinahan, 2016). While extensively studied in oncology, radiomics has been underexplored in the context of carotid plaque analysis (DeMarco and Porambo, 2020; Saba et al., 2024; Vacca et al., 2024). By extracting texture features from MR images, machine learning classifiers can be developed to predict carotid plaque vulnerability, potentially paving the way toward personalized treatment strategies (Le et al., 2021; Meddings et al., 2024; Lambin et al., 2017).

While radiomic image texture features may help characterise lesion phenotypes, further insight into plaque stability may be provided by the biomechanical environment (Curcio et al, 2023; Mahmoud, 2019; Oyejide, 2023). Hemodynamic forces – particularly wall shear stress (WSS) - have long been associated with the destabilization of established plaques near regions of complex blood flow, such as the carotid bifurcation and regions of stenosis (Duerinckx and Shaaban, 2000). However, the exact role of WSS in both early-stage and late-stage atherosclerosis is largely unresolved (Gijsen et al, 2019), and it remains inconclusive as to whether low and/or oscillatory WSS triggers elevated, atherosclerosis-exacerbating, low-density lipoprotein (LDL) levels (Peiffer et al, 2013; Weinberg, 2022).

Another important aspect of the mechanics environment is the vessel structural stress (VSS), which represents the intra-mural forces acting within the plaque. These forces are associated with plaque rupture, which occurs when VSS overcomes the material strength of the arterial wall. This can lead to a stroke or TIA via an exposed necrotic core and subsequent thrombosis (Bentzon, 2014). Therefore, regions of high VSS may be associated with plaque vulnerability due to the potential for lesion rupture (Sadat et al, 2010; Tang et al, 2009).

This study aims to integrate vessel wall mechanics with radiomics-based texture analysis into a unified mechano-radiomic model for predicting plaque vulnerability. Mechanical stress maps derived from 3D computational models derived from multi-slice MR images are analysed alongside radiomics features extracted from the original imaging data. This combined approach seeks to enhance the understanding and prediction of plaque vulnerability, ultimately contributing to improved risk stratification and patient management strategies.

2. Materials and Methods

2.1. Patients and High-Resolution, Multi-Contrast MRI

Data from sixty-five patients from the CARE-II multicentre clinical trial were used in this study (ClinicalTrials.gov Identifier: NCT02017756) (Zhao et al, 2017). Patients in this study cohort had symptoms attributed to one of their carotid arteries, such that the dataset could be split evenly into culprit and non-culprit lesions. Culprit lesions were defined as the side deemed responsible, following clinical assessment and image analysis, for a transient ischaemic attack (TIA) or stroke. A summary of cohort demographics is provided by Table 1.

Patients underwent bilateral MRI of their carotid arteries with the following sequences: three-dimensional (3D) time-of-flight (TOF), T1-weighted (T1w) quadruple inversion recovery (QIR), T2-weighted (T2w) double inversion recovery (DIR), and Magnetisation Prepared - Rapid Acquisition Gradient Echo (MP-RAGE). Detailed MR parameters can be found in the reference paper for the CARE-II trial (Zhao et al, 2017).

2.2. CFD and FEA Methodology

The mechanical environment of the artery was simulated using a one-way fluid-structure-interaction (FSI) approach within the ADINA (Bentley Systems Incorporated, Exton, PA) software. Computational fluid dynamics (CFD) simulations were used to estimate wall shear stresses (WSS) and endothelium pressures. Pressure results were then carried forward from the CFD simulation to a finite element model ahead of a structural deformation finite element analysis (FEA) simulation. Forces derived from fluid dynamics were applied as boundary conditions in the structural model, without feedback from the structural simulation result.

The models which were used for CFD and solid deformation analysis were extrapolated from manually segmented contours of the lumen, outer wall boundaries and various atherosclerotic components including intra-plaque haemorrhage (IPH), lipid core and calcification. Segmentations were performed manually by with reference to T1-weighted, T2-weighted, TOF, and magnetisation-prepared rapid gradient-echo (MP-RAGE) image modalities. Segmentations were performed on a minimum of 16 axial slices and verified by a minimum of 2 reviewers. This process was supported by a custom-designed software (CASCADE; University of Washington, Seattle, USA) (Zhao et al, 2017).

Segmentations were conducted in a custom-built MATLAB interface (The MathWorks, Natick, MA). The endothelium surface was modelled by discrete 2D shell elements for CFD analysis. Using the finite element pre-processing software Hypermesh (Altair Engineering, Inc., USA), given the complexity of the geometries, 3D models were generated with an unstructured mesh protocol, consisting of tetrahedral elements. Proximity-based mesh density adjustments were made to allow for a greater number of elements in tight regions, such as thin walls. While 3D meshes were bespoke for each patient and artery side, the total number of elements were of the order 10⁵. Grid-independence validation tests were done by adjusting the number of elements and checking for wall pressure convergence from identical CFD simulation conditions.

The incompressible Navier-Stokes equations were solved using a Lagrangian-Eulerian formulation in ADINA to simulate blood flow, where $\rho$ represents blood density, $u$ velocity, $t$ time, $p$ pressure and $u$ viscosity.

$$\rho {\displaystyle \frac{\partial u}{\partial t}}+\rho \left[u.\nabla u\right]=-\nabla p+{\mu }^{2}u$$

CFD simulations were performed based on time-varying blood pressure inlet and outlet boundary conditions, to simulate a cardiac cycle and an estimation of time-varying wall shear stress (WSS). WSS was computed at discrete points on the inner wall. The results were subsequently transposed onto a 2D mechanics image through a flattening procedure to provide a spatial representation. Blood pressure data was obtained from patient-specific systolic and diastolic measurements from the CARE-II dataset.

Wall pressures from the CFD simulations were applied to compute the intramural stresses and deformations in the artery. Slice-wise vessel structural stress (VSS) was calculated from finite element analysis (FEA). Stress measurements were based on local maxima observed during the cardiac cycle, capturing the peak structural stress experienced by the vessel.

The principal stress in the first direction (Stress-P1) was used to represent VSS, as it corresponds to the maximum tensile stress invariant. This stress is aligned with the material fiber orientation, which is indicative of fiber material strength and structural integrity. To model the arterial wall's behaviour accurately, a modified Mooney-Rivlin strain energy density function (SEDF) was employed. This provided a detailed distribution of structural stress across the artery, accounting for its non-linear elastic properties.

2.3. Boundary Conditions, Parameters and Assumptions

The computational fluid dynamics (CFD) and finite element analysis (FEA) simulations were based on several key boundary conditions, parameters, and assumptions:

• Inlet and outlet boundary conditions: Time-varying inlet pressure profiles based on patient-specific measurements were applied at the proximal end of the arterial segment. This pressure ranged between the patient’s diastolic and systolic blood pressure measurements, which were provided as part of the CARE-II dataset. A relatively small distal resistance assumption justified a grounded outlet pressure boundary condition.
• Blood flow characteristics: Blood flow was modelled as incompressible, viscous, and laminar. The blood was assumed to behave according to Newtonian fluid dynamics, which is a reasonable assumption for larger arteries. The viscosity of blood was set to 3.5 × 10−3 Pa·s, and the density was defined as 1000 kg/m³, reflecting typical physiological values for human blood.
• Wall conditions: The wall was modelled as rigid, and a no-slip boundary condition was imposed at the arterial walls. This means that the velocity of the blood at the wall was set to zero. The inner wall experienced loading defined by pressures extracted from the CFD simulations. Elements were free to move in radial and circumferential directions during loading.
• Solver and convergence settings: Transient CFD simulations were conducted under time-varying inlet and outlet pressures to simulate a cardiac cycle. An implicit Navier-Stokes solver was employed using the ADINA-F module. For FEA, an implicit, full-Newton, large strain and large displacement formulation was used within the ADINA-S module. Convergence criteria were set to 10^-6 for both energy and displacement residuals. Time-steps of 0.009s for 0.9s long loading cycles were chosen following time-step independence testing to balance computational time with precision. Grid-independence testing was also employed to determine a mesh size which addressed a similar computational trade-off.
• Material properties: The carotid artery material’s hyper-elastic constitutive model which relates stress to strain was governed by a modified Mooney-Rivlin strain-energy density function (SEDF). The Mooney-Rivlin model has consistently been shown to provide adequate approximation to the behaviour of the artery wall when compared with other constitutive models (Teng et al, 2015). Finite element analysis (FEA) was employed to solve for stress. The SEDF was modified for each type of biological material: wall; lipid-rich necrotic core (LRNC); intra-plaque haemorrhage (IPH); fibrous tissue; and calcification. Properties were chosen based on previous uni-axial tension experiments (Teng et al, 2014).

These boundary conditions and assumptions were selected to balance computational efficiency with physiological realism, ensuring that the results provide insight into both fluid dynamics and structural stress without overly complicating the model.

2.4. Mechanics Post-Processing

Centrelines were extracted from the 3D artery models using the vascular modelling toolkit VMTK (Antiga et al, 2008). An inward remodelling (IR) index was calculated at each artery wall node by dividing its distance to the centreline by a reference length in the common carotid artery (CCA). The single stenosis measurement, typically the ratio of local diameter to reference diameter, is related but not equivalent to a local maxima on the inward remodelling map. Similarly, the wall shear stress (WSS) was calculated at each node on the inner wall. WSS and IR index sampling at each node was non-uniform due to uneven node clustering. To address this, a patching procedure was developed to convert the wall into a uniform set of discrete zones that included the common and internal carotid arteries. Each discrete zone reported the average values from the nodes within it.These zones can be interpreted as the pixels of a 2D map generated from the mechanics data. A radiomics-style analysis was then applied to these maps to evaluate the predictive potential of the IR index and WSS maps.

Vessel structural stress was calculated at critical regions of the artery wall, where localized high stress concentrations occurred, such as the bifurcation and plaque shoulder. VSS feature extraction focused on the single value of the critical VSS. The cross-section selected for analysis was the one where the critical VSS occurred—typically at local maxima at the bifurcation or plaque shoulder region. The patching and mapping procedure is further outlined in Figure 1, which presents a WSS map used for mechanics map feature extraction. Inward remodelling (IR) index maps were based on the same principle. A slice-wise VSS distribution, approximately perpendicular to the lumen centerline, is illustrated in Figure 2.

2.5. Radiomic Feature Extraction

The open source PyRadiomics package (v2.2.0) (van Griethuysen et al, 2017) was used to extract averaged radiomic features from multiple MR image slices, as well as from a single 2D mechanics map, for each artery. Features belonging to the following groups were determined: first-order intensity histogram statistics; Gray Level Co-occurrence Matrix features (GLCM) (Dinstein, Haralick and Shanmugan, 1973; Haralick, 1979; Rundo et al, 2019); Gray Level Run Length Matrix (GLRLM) (Galloway, 1975); Gray Level Size Zone Matrix (GLSZM) (Angulo, Meyer and Thibault, 2014); Gray Level Dependence Matrix (GLDM) (Sun and Wee, 1982); and Neighbouring Gray Tone Difference Matrix (NGTDM) (Amadasun, 1989). The default PyRadiomics values of 25 for fixed bin width, and 107 total features extracted, were used. No resampling was applied.

This study used features based on the T2-weighted modality. It has been shown (Meddings, 2024) that T2-weighted radiomic features may offer better sensitivity to variations in the region of interest.

2.6. Predictive Modelling

Logistic regression models were developed to classify arteries into culprit and non-culprit categories, using a combination of radiomics features and biomechanical parameters. The modeling process was designed to ensure rigor and reproducibility, employing both feature selection and robust cross-validation techniques.

Prior to modeling, features underwent pre-processing to ensure consistency and robustness. Radiomic features with poor intra-class correlation coefficients (Meddings, 2024) were excluded to minimize noise and variability. Remaining features were standardized to have zero mean and unit variance to prevent dominance by features with larger scales.

Feature selection was conducted using the least absolute shrinkage and selection operator (LASSO) regression. This analysis was used to develop sparser models and reduce over-fitting (Tibshirani, 1996), allowing for the retention only of those features with strong predictive capabilities. Hyperparameters for LASSO were optimized through cross-validation within the training folds to maximize predictive accuracy while maintaining model interpretability.

Combining five-fold cross validation and leave-out-out cross-validation (LOOCV) (Cheng, Fernando and Garrick, 2017) provided a comprehensive perspective on model performance. The application of LOOCV ensured that each test case remained unseen during model training without leakage, while the size of the training dataset was maximised. Five-fold cross-validation splits offered insights into the variability and generalizability of the models. These complementary validation strategies bolster confidence in the robustness of the mechano-radiomic model.

LOOCV was used to evaluate model performance, while five-fold cross-validation was used to generate receiver-operating characteristic (ROC) curves, and the area under the curve (AUC) for ROC evaluation. Model weights presented in subsequent analysis were average across all folds in the LOOCV scheme, which benefited from more training data – an important factor given the small cohort size. The overall schematic is illustrated in Figure 3 below.

2.7. Statistical Analysis

Results were considered statistically significant if p < 0.05. McNemar’s test for paired samples was used to compare error rates between the two classifiers. The null hypothesis is that the error rates are equal (Dietterich, 1998). Significance between two ROC curves were assessed using the DeLong test (DeLong, 1988).

3. Results

3.1. Model Evaluation

The evaluation of the predictive models was performed using both leave-one-out cross-validation (LOOCV) and five-fold cross-validation, as summarised in Table 2 and Table 3. In Figure 4, the receiver-operating characteristic (ROC) curves are plotted for the radiomics, mechanics, and combined mechano-radiomic models. These curves illustrate the trade-off between sensitivity and specificity for each model.

In summary, both mechanics [accuracy = 0.68, AUC = 0.75 $\pm$ 0.03] and radiomics models [accuracy = 0.67, AUC = 0.71 $\pm$ 0.04] showed greater predictive capabilities for culprit lesions than the conventional measurement of luminal stenosis alone [accuracy = 0.62, AUC = 0.57 $\pm$ 0.04] (p < 0.001)]. The combined mechano-radiomic model [accuracy = 0.76, AUC = 0.82 $\pm$ 0.04] provided a significant improvement over the radiomic-only model (p = 0.037).

3.2. Feature Selection

Models were built by minimising over-fitting using LASSO feature selection. This approach penalised less informative features, yielding sparser models can be more easily interpreted. The co-efficients of the more predictive features, as well as the results of some less predictive but commonly cited texture features from the radiomic, mechanics and mechano-radiomic models, are illustrated in Table 4, Table 5 and Table 6 respectively, based on the LOOCV scheme. Feature weights are presented as averages across all folds, along with a 95% confidence interval.

A mechano-radiomic signature may be represented by the top 10 weighted features across all cross-validation folds during model evaluation, which combined features from mechanics and radiomics models. This is presented in Table 7 as a set of features and odds ratios, derived from the logistic regression model weights $\beta .$ The odds ratio, defined as $e^{\beta }$, represents the predicted change in the odds of the outcome for a unit increase in the normalized predictor variable, if all other variables were to be held constant.

4. Discussion

The performance of models combining traditional medical image-based radiomic features with features derived from mechanics maps emphasises that integrating radiomics and a comprehensive biomechanical environment may significantly improve identification of culprit lesions compared with more conventional methods.

The predictive model contributions of shear stress (WSS) are particularly noteworthy: while maximum time-averaged wall shear stress (WSS) appeared to be a predictor of non-culprit plaques, the 3D distribution of the WSS map (WSS Zone Percentage) and the inward remodelling (IR) interquartile range measurements did demonstrate strong predictive capabilities, supporting the notion that distribution metrics of mechanical properties may serve as stronger predictors than single-point measurements, such as local minima or maxima. Although WSS might be more relevant to atherosclerosis development in earlier stages, it may still play a crucial role in determining plaque vulnerability in advanced stages of the disease.

Critical vessel structural stress (VSS) did not exhibit a strong multivariate odds ratio, as indicated by its low logistic regression coefficient. This is likely due to its dependence on patient-specific blood pressure levels. The culprit vessel may not necessarily experience higher stress for several reasons: tissue more fatigued under cyclical loading; more extensive vascular remodelling; ipsilateral plaque material strength is weaker; tissue is more erosion-prone; tissue has concentrations of micro-stresses; and other factors unaccounted for by the simulations.

The radiomic model in Table 4 is composed mainly of features which relate to image complexity, such as Local Homogeneity (GLCM), Small Dependence Low Gray Level Emphasis (GLDM), and the Local Homogeneity Idmn (GLCM) feature which was featured in the combined model signature in Table 7. The combined model utilises a wide range of texture features across the different image modalities, balancing size zone features of the WSS maps such as zone percentage (number of distinct WSS zones) with inward remodelling (IR) index map texture features such as Dependence Variance (GLCM). As with individual models, higher-order image texture features consistently had high weightings, indicating the predictive capabilities of radiomics in general. Our reported accuracy and AUC for radiomics are in line with similar studies. Using (Le, 2021) as an example, radiomics predictive capabilities on CT carotid images had mean AUC across 5 folds of 0.67, while the corresponding value in our study was 0.73. The improvement when combining radiomics with other biomarkers in this paper is demonstrated by an AUC of 0.73, so the combined AUC measurement of 0.82 is a positive development towards having a more comprehensive model of carotid risk management using AI-assisted high-resolution medical imaging.

The study is limited by its relatively small patient cohort and its retrospective design. The clinical goal of this research is to enable proactive diagnosis before a cardiovascular event occurs, with the potential to integrate it into screening protocols for both symptomatic and asymptomatic patients. An implicit assumption underlying this work is that the post-event mechanical and imaging structures resemble those present before the event.

Future work will be directed towards performing similar analyses for larger patient numbers, with an external test set to further validate the mechano-radiomic model. For example, selected features would benefit from stability analysis across larger cohorts or different populations. The presented model should be combined with further biomarkers, such as from functional MRI imaging, and deeper geometrical descriptions of arteries. Eventually, a deep learning method based on convolutional neural networks may facilitate a more efficient transition from imaging to mechanical analysis.

5. Conclusions

In conclusion, the mechano-radiomic model significantly improved plaque risk assessment over either of its constituent mechanics and radiomics methods alone, as well as the conventional stenosis measurement. Thus, in the context of carotid artery disease, combining image texture analysis from MR images with the mechanical environment offers a novel approach for risk stratification in plaque vulnerability assessment.