The study included 443 patients from two academic centers with a confirmed diagnosis of AIS, due to large vessel occlusion in the anterior circulation [internal carotid or middle cerebral artery (MCA)] confirmed on computed tomography angiography (CTA) who underwent mechanical thrombectomy between 2014 and 2020. All included patients underwent standardized acute stroke imaging that includes non-contrast head CT and CTA of the head and neck. Patients without CTA or with any missing pertinent clinical or radiological data were excluded. The primary clinical outcome of interest was the modified Rankin score at 90 days (mRS-90). The study protocol was approved by the Institutional Review Board. The data that support the findings of this study could be available from the corresponding author upon reasonable request.

Collected clinical and demographic characteristics were age, sex, baseline National Institutes of Health (NIH) Stroke Scale (NIHSS), time from symptom onset to admission, and comorbidities (diabetes, hypertension, hyperlipidemia, previous stroke, and cardiovascular disease, such as myocardial infarction and arrhythmia), in addition to blood glucose and blood pressure levels. Interventional and post-interventional features, such as modified treatment in cerebral infarction (mTICI) score informing reperfusion status, were excluded from the models in line with the main purpose of the study.

Quantitative imaging feature extraction was performed using e-Stroke software (Brainomix, Oxford, UK) for automated calculation of Alberta Stroke Program Early CT Score (ASPECTS; e-ASPECTS) and estimated acute infarct volumes on non-contrast CT (14–17). e-ASPECTS uses ML classification to distinguish and segment regions that contain signs consistent with the acute ischemic change in order to output both (total and per ASPECTS region) and total e-ASPECTS volumes. Additional novel features that were extracted using e-Stroke software included non-acute infarct volume, total brain volume, and atrophy, which were quantified using cerebrospinal fluid (CSF) segmentation volumes in both the lateral ventricles and the cortical sulci separately and expressed as percentages. e-CTA (Brainomix, Oxford, UK) identifies large vessel occlusion site and quantifies the volume of collateral circulation deficit both as a percentage of the total volume and using the CTA collateral score (CTA-CS) (18–20). Novel experimental outputs from e-CTA included the absolute volume of the vessel density deficit in MCA territory relative to the contralateral hemisphere.

Baseline features were categorized into clinical and imaging feature groups. Standardization scaling of continuous and ordinal feature values was applied to obtain a mean of zero and a standard deviation (SD) of 1, in order to facilitate the algorithm learning process and improve the prediction results. Random splitting of the datasets into a training set and a testing set was applied with a ratio of 75–25%, respectively. The mRS-90 was dichotomized with mRS 0–2 representing a good functional outcome.

The features were divided into four subsets, which are as follows: (1) clinical features, (2) imaging features, (3) combined clinical and imaging features, and (4) selected features. A model-based approach was applied using sequential backward feature selection with a bagging classifier, where an algorithm sequentially removes features from the full feature set until the removal of further features decreases the classifier performance.

Statistical assessment of each clinical and image-based feature in relation to mRS-90 was assessed using the chi-square test for the categorical variables, Wilcoxon Rank-Sum test, and t-test for the ordinal and continuous variables depending on the normality of their distributions. Statistical analysis was done using Python (version 3.9) and the SciPy library. Values of p < 0.05 were considered statistically significant.

For the purpose of mRS-90 prediction, supervised ML classification methods were deployed. The ML algorithms used were as follows: k-nearest neighbors, random forests (RF), gradient boosting (GB), and Extreme Gradient Boosting (XGBoost). The models were constructed using the Scikit-learn library.

As a first step, k-fold cross-validation of 10-folds was performed during the training for each model, which divides the training set into 10 subsets (9 for training and 1 for validation), where the training and validation sets change and iterate over the 10-folds. The model hyperparameters were optimized by means of a grid search approach, where for every model and for each hyperparameter a set of possible values was manually defined and evaluated exhaustively in every iteration to determine the values corresponding to the model's highest performance, with an area under the receiver operating characteristic curve (AUC) as scoring metric. The ML models were trained using each of the 4 different feature categories. Subsequently, for every feature group, we tested the models' performance on the testing set of patients.

From the output of the grid search, the best performing model was chosen. Finally, automatic Bayesian hyperparameter tuning with the Optuna framework was used on the best performing model to boost its performance and achieve finer tuning. The evaluation metrics used were accuracy, F1 score (for mRS-90 ≤ 2 and mRS-90 > 2 predictions), and AUC.

To enhance the model's explainability and perform a feature importance analysis, we used the method SHAP, which is based on game theory and consists of computing Shapley values reflecting the contribution of each feature in the predictions of the model (21–23). The method allows for the identification of features with the most influence on model output and measures the impact if each variable was to be removed while taking into account the interaction with other variables that provide insight on the relative importance of the features used by the model for its prediction decision process.

An illustrative summary of the methods is provided in Figure 1.

Of 443 total patients (266 from the first center and 177 from the second center), 293 patients met the study inclusion criteria. The remainder of the patients were excluded for lacking relevant variables mainly CTA images. In total, 101 patients had a favorable functional outcome (mRS-90 ≤ 2), while 192 patients had unfavorable functional outcomes (mRS-90 > 2). The median age of included patients was 71 years and 49% (n = 143) were women (Table 1).

Favorable clinical outcome was significantly associated with younger age (p < 0.0001), female sex (p = 0.016), and lower baseline NIHSS score (p < 0.0001). Patient comorbidities were not significantly different between outcome groups (Supplementary Table S1) and therefore were not included in the development of ML models.

Non-contrast CT imaging features associated with favorable outcomes included greater e-ASPECTS (p = 0.002), larger brain volume (p = 0.043), smaller cortical CSF volume (p < 0.0001), smaller lateral ventricle volume (p < 0.0001), smaller acute ischemic volume (p = 0.047), and non-acute ischemic volume (p = 0.040). Collateral circulation deficit volume on CTA was significantly lower in the favorable outcome group (p = 0.001; Table 1).

Performance evaluation of each ML model following grid search optimization is presented in Supplementary Table S2. The selected features were as follows: baseline NIHSS, age, occlusion side, local M5 infarct volume, local lentiform infarct volume, brain volume, percentage of lateral ventricle volume, collateral vessel deficit volume, and the time interval from symptoms onset to admission.

We calculated accuracy and AUC for each of the feature groups on the training data set. For clinical data, the model with the highest AUC score was the XGBoost classifier (XGB) with an AUC of 81%. For imaging features, the best model was XGB at 79% AUC. For combined clinical and imaging features, the best model was also XGB with an AUC of 80%. Selected features yielded an AUC of 84% (Figures 2A–D).

The XGB model was selected for further optimization as it consistently achieved a high performance in the four feature groups and had the highest overall AUC scores. Using Bayesian hyperparameter tuning with a stratified cross-validation of 10-folds to refine the XGB model, the final performance metrics on the testing set of patients were AUC = 84%, accuracy = 77%, F1-score (mRS ≤ 2) = 67%, and F1-score (mRS > 2) = 82% for the selected features. The final results are shown in Figures 2E,F.

Following prediction modeling of mRS-90, a feature importance rank for the patient cohort was established by calculating SHAP values for XGBoost which revealed that the top indicators of clinical outcome prediction for the model were by order of importance, i.e., age, baseline NIHSS, occlusion side, cortical CSF volume percentage, lateral ventricle volume percentage, e-ASPECTS, and circulation deficit volume (Figure 3E). The overall impact of each feature is represented by feature SHAP general values as shown in Figure 3F and model predictions were able to be reviewed and assessed regarding each predictor for each patient instance. The SHAP force plot allows for an interactive visualization of all the study populations clustered by their feature value similarity and ranging according to their specific model output (Figures 3A,B). Individual patient predictions can be extracted to visualize which features played a role in their classification and what their feature values were. Examples of predictions for a patient with poor outcomes and a patient with favorable outcomes are shown in Figures 3C,D, respectively.

Limitations to our study may relate to the population size, which is although relatively large as a two-center study, similar to most ML studies, it could benefit from a larger cohort size for ML purposes. Datasets with diverse origins and a higher number of participants are warranted to further validate the robustness of the models for future generalizability on independent cohorts. In addition, with the absence of consistent information on pre-morbid functional status, we have not included this variable. Future planned steps exist for validating these tools prospectively and on larger multi-center datasets for further optimization of this approach.