INTRODUCTION
Ossification of the posterior longitudinal ligament (OPLL) is a degenerative spinal condition characterized by pathological ossification along the posterior longitudinal ligament, which can lead to significant spinal canal stenosis, neurological deficits, and, in severe cases, cervical myelopathy [1]. OPLL most commonly affects the cervical spine, with a cross-sectional study reporting a prevalence of approximately 4.1% in the Chinese population [2]. The progression and clinical impact of OPLL can vary significantly among patients, creating substantial challenges in management and prognosis.
The natural history of OPLL is not fully understood, but it is known that the ossified mass can enlarge over time, exacerbating spinal cord compression and worsening clinical outcomes [3]. Some studies have shown that more than 60% of patients experience OPLL progression after posterior cervical surgery [4-7]. This progression can negatively impact the outcomes of spinal decompression surgery, leading to symptom recurrence or the need for additional surgical interventions [8]. Therefore, predicting the likelihood of postoperative progression is crucial for optimizing treatment plans and improving long-term patient outcomes. While some clinical studies have identified risk factors for postoperative OPLL progression, such as age, OPLL subtype, and surgical technique, no studies have yet developed predictive models for this progression [9-12].
Radiomics, an emerging field that extracts a large number of quantitative features from medical imaging, has shown promise in enhancing the prediction of various disease outcomes [13]. By analyzing these high-dimensional data, radiomics can capture intralesional heterogeneity that might not be apparent on standard imaging. In recent years, radiomics has been increasingly applied in oncology for tumor characterization, treatment response prediction, and survival analysis, but its application in degenerative spinal diseases remains underexplored [14-17]. Some researchers have applied this approach to predict the minimal clinically important difference after OPLL surgery, demonstrating the potential of radiomics in OPLL prognostication by providing additional information for predicting surgical outcomes through the extraction of high-dimensional features from OPLL lesions [18].
This study aims to develop interpretable radiomics models to predict the postoperative progression of cervical OPLL. By integrating clinical and radiomics data, the study seeks to enhance the accuracy of prognostic predictions, potentially leading to more personalized treatment and postoperative management strategies for patients undergoing spinal surgery for OPLL.
MATERIALS AND METHODS
1. Patients
The study was performed according to the Helsinki Declaration (http://www.wma.net/en/30publications/10policies/b3/) and approved by the Institutional Review Board (IRB) of Peking University Third Hospital on June 18, 2024 (approval No. IRB00006761-M2024489). Patient informed consent was waived. Due to the retrospective study design and the use of existing data, the requirement for informed consent was waived. Patient information was retrieved from electronic medical records for a retrospective analysis of individuals diagnosed with OPLL at Peking University Third Hospital between October 2006 and September 2022. Inclusion criteria included: (1) patients who underwent posterior spinal surgery; (2) availability of at least 2 computed tomography (CT) examinations spaced at least 1 year apart. Exclusion criteria encompassed: (1) patients who received additional anterior procedures or revision surgeries and (2) Initial OPLL lesion has too few voxels for feature extraction. CT scans were performed both before surgery and at the final follow-up. Demographic data (such as sex and age) were extracted from the electronic health records, and radiological assessments were made on the x-ray, CT, and magnetic resonance images from both the preoperative and the final postoperative visits. The patient inclusion and exclusion process are illustrated in Fig. 1.
2. CT Acquisition Scan Protocol
Initial cervical spine CT scans were performed using the following equipment: GE CT620 (GE Healthcare, Chicago, IL, USA), GE CT750HD (GE Healthcar), GE Rev (GE Healthcare), SE FLASH (Siemens Healthineers, Erlangen, Germany), UI CT790 (United Imaging Healthcare, Shanghai, China), UI CT960 (United Imaging Healthcare), UI CT 710 (United Imaging Healthcare). All CT scans were performed using parameters set at 120 kV and smart mA, with reconstructions oriented along the long axis of the cervical pedicles. The slice thickness was set at 3 mm.
3. Calculation of the Annualized Growth Rate of OPLL Volume
First, the preoperative and the most recent postoperative follow-up axial CT images of the cervical spine were downloaded from the hospital’s PACS (Picture Archiving and Communication System). These images were then imported into ITK SNAP (ver. 4.0.1) software. The window width and level were adjusted to optimally display the lesions. An experienced radiologist with 4 years of experience segmented the OPLL lesions at the C1–7 levels layer by layer on the images, marking them as the region of interest (ROI). After the initial segmentation, a senior musculoskeletal radiologist with 22 years of experience reviewed and modified the ROI. Upon completion of the segmentation, the software automatically generated a 3-dimensional (3D) volume of interest (VOI) and calculated the lesion volume. The annual growth rate of OPLL volume was calculated using the following formula:
According to Katsumi et al. [11], an annual growth rate greater than 7.5% was defined as the OPLL progression (P) group, and an annual growth rate less than 7.5% was defined as the nonprogression (NP) group. Fig. 2. illustrates the segmentation of OPLL and the calculation of the annualized growth rate.
4. Radiomics Feature Extraction
Feature extraction was performed on the VOI of the preoperative OPLL volume. Due to the use of multiple CT devices in this study and to minimize the impact of imaging protocol variations, the voxels of the original images were resampled to a resolution of 1 mm³× 1 mm³× 1 mm³ prior to feature extraction. The min-max normalization method was then used to normalize the pixel values of the images to the [0, 1] range. CT-based radiomic features were extracted using the Pyradiomics package (ver. 3.0) in Python (ver. 3.7.3), following the guidelines of the Image Biomarker Standardization Initiative [19].
The manually crafted features can be divided into 3 groups: (1) shape (14 features), (2) intensity (18 features), and (3) texture (73 features). Shape features describe the 3D shape characteristics of OPLL. Intensity features describe the first-order statistical distribution of voxel intensities within the OPLL. Texture features describe patterns or the second- and higher-order spatial distributions of intensities. Various methods were used to extract texture features, including the Gray Level Co-occurrence Matrix, the Gray Level Run Length Matrix, the Gray Level Size Zone Matrix, and the Neighboring Gray Tone Difference Matrix methods.
In addition to extracting features from the original images, the same types of features were extracted from filtered images using several filters, including the Laplacian of Gaussian (with sigma values of 1.0, 2.0, and 3.0 mm), Wavelet, Local Binary Patterns in 3D, Exponential, Square, SquareRoot, Logarithm, and Gradient filters. In total, 1,834 radiomic features were extracted from each VOI.
5. Feature Selection
In the feature selection process, we first group the dataset and split it into training and testing sets in a 7:3 ratio. Following this, we standardize the features to have a mean of zero and a standard deviation of one, ensuring that all features contribute equally to the analysis and preventing larger scale features from disproportionately influencing the model. The standardized dataset then undergoes an initial feature removal step where constant features are discarded. Next, we address multicollinearity by calculating the correlation matrix for all numerical features. We identify pairs with correlation coefficients greater than 0.8, retaining only one feature from each pair and discarding the others to remove highly correlated features. Subsequently, we applied the least absolute shrinkage and selection operator (LASSO) for further feature selection, using 5-fold cross-validation to determine the optimal regularization parameter (alpha). To further reduce feature dimensionality and minimize interfeature correlation, we employed principal component analysis (PCA) to extract principal components that cumulatively explained 95% of the variance, enhancing the model’s computational efficiency and robustness.
6. Collection of Baseline Clinical Information
Age, sex, surgical approach, number of treatment segments, initial OPLL volume, and OPLL type were recorded as baseline clinical and radiological information. The surgical approaches included laminoplasty (LMP) and posterior decompression with instrumented fusion (PDF). OPLL types included continuous, segmental, mixed, and localized types. On cervical spine x-rays, the C2–7 angle and range of motion (ROM) are measured. The C2–7 angle is defined as the angle between the lower endplates of C2 and C7. Extension ROM is calculated as the change in the C2–7 angle from the neutral position to the extension position, while flexion ROM is the change in angle from the neutral to the flexion position. On sagittal MRI T2-weighted images, the signal intensity at the narrowest segment of the spinal cord is classified into 3 grades: grade 0: no high signal intensity; grade 1: mild high signal intensity (blurred signal); grade 2: high-intensity signal (bright signal).
7. Statistical Analysis
1) Statistical methods for clinical variables
Continuous variables, including age, initial OPLL volume, C2–7 angle, extension ROM, and flexion ROM, were tested for normality (Shapiro test) and subsequently analyzed using either the t-test or the Mann-Whitney U-test. Categorical variables, including sex, surgical approach, number of treatment segments, OPLL type, and intramedullary signal intensity (ISI), were analyzed using chi-square tests. We employed the Kruskal-Wallis H-test to evaluate the significant differences in annual progression rates among different OPLL types. Following this, Dunn test was conducted to identify specific pairs of OPLL types with significant differences. This analysis utilized pairwise Mann-Whitney U-tests with Bonferroni correction to control for type I errors resulting from multiple comparisons.
2) Development of logistic regression models
Spearman correlation analysis was used to calculate the correlation coefficients between different principal components, which were visualized using a heatmap. The selected principal components were used to develop a logistic regression (LR) model, and the optimal model parameters were determined through grid search. The Rad-score was calculated based on the principal components and the coefficients of the LR model, as described in Supplementary Material 1. Clinical variables with p<0.05 in univariable analysis were used to develop a clinical model. A combined model was constructed using both clinical variables and the Rad-score. The model’s performance was evaluated using the area under the receiver operating characteristic curve (AUC) with a 95% confidence interval (CI), and the optimal classification threshold was determined using the maximum Youden index. Sensitivity, specificity, precision, recall, F1 score, accuracy, positive predictive value (PPV), and negative predictive value (NPV) were calculated at this threshold.
3) Model evaluation, clinical utility, and interpretability techniques
A nomogram was constructed based on the LR model to provide an individualized prediction tool for clinical use. Calibration curves were used to assess the agreement between predicted probabilities and observed outcomes for the LR models. Decision curve analysis (DCA) was conducted to determine the clinical utility of the model by quantifying the net benefit across a range of threshold probabilities. SHAP (SHapley Additive exPlanations) analysis was performed to interpret the contributions of individual features to the model’s predictions. SHAP values indicate the impact of each feature on the prediction, allowing for a comprehensive assessment of feature importance and how different features influence the model output.
A 2-sided p-value of less than 0.05 was considered statistically significant.
The radiomics workflow in this study is shown on Fig. 3.
RESULTS
1. Baseline Patient Information
A total of 473 patients were included in this study, with a follow-up period ranging from 12 to 202 months and a median follow-up time of 22 months. Among these, 191 patients experienced OPLL progression and were classified into the P group (annual progression rate > 7.5%), while 282 patients did not show OPLL progression and were classified into the NP group. The P group had a significantly higher average initial OPLL volume (2250 mm³) compared to the NP group (1,883 mm³, p=0.027). Additionally, the P group demonstrated a significantly higher median annual progression rate of 13.7%, while the NP group had a median annual progression rate of only 2.0% (p<0.001). The P group had a significantly lower average age (52.94±8.89 years) compared to the NP group (57.06±9.00 years, p<0.001). For preoperative OPLL type, the distribution differed significantly between the 2 groups (p<0.001). In terms of other preoperative clinical and radiological information, there were no significant differences between the 2 groups in sex, C2–7 angle, extension ROM, flextion ROM, ISI, number of treatment segments, surgical approach, and preoperative OPLL volume. Detailed patient information is shown in Table 1.
2. Analysis of Progression Rates in Different Types of OPLL
In this study, the mixed type of OPLL exhibited the highest median annual progression rate at 7.5% (interquartile range [IQR], 4.3%–14.9%), followed by the continuous type at 5.1% (IQR, 1.8%–9.8%) and the segmental type at 4.9% (IQR, 0%–12.5%). The localized type showed the lowest median annual progression rate of 0% (IQR, 0%–4.3%) (Fig. 4A). The Kruskal-Wallis H-test revealed significant differences in annual progression rates among the different OPLL types (H=39.703, p<0.001). Post hoc pairwise comparisons using the Bonferroni correction demonstrated that the mixed type had significantly higher progression rates compared to the segmental (p=0.002), localized (p<0.001), and continuous types (p=0.046). Additionally, the progression rate of the localized type was also significantly lower compared to the segmental type (p=0.001) and the continuous type (p<0.001) (Supplementary Tables 1, 2). These findings underscore the variability in progression rates across OPLL types, suggesting that specific type-related factors may influence disease progression.
During the follow-up period, regarding OPLL classification, 5 patients progressed from localized type to segmental or mixed type, 42 patients progressed from segmental type to mixed or continuous type, 30 patients progressed from mixed type to continuous type, and 1 patient with continuous type developed a new segment, becoming mixed type (Fig. 4B).
3. Feature Selection and Model Development
Following correlation coefficient analysis and LASSO, 16 radiomics features were retained (Supplementary Table 3). After PCA, the first 13 principal components were retained. Detailed procedures and parameters for LASSO and PCA are provided in Supplementary Materials 2.1 and 2.2. As shown in Supplementary Fig. 4, the correlations between the 13 principal components were very low, making them suitable for LR modeling. The calculation formula for the Rad-score is provided in Supplementary Material 1. The clinical model was developed using 3 variables: age, postoperative OPLL type, and preoperative OPLL volume. The combined model was constructed using both the Rad-score and clinical variables. Table 2 summarizes the performance metrics of the Rad-score, clinical, and combined models on the training and testing sets. The combined model achieved the best overall performance, with the highest AUC for both the training set (0.777; 95% CI, 0.729–0.826) and the testing set (0.751; 95% CI, 0.670–0.828). The AUC of the combined model on the testing set was significantly higher than that of the Rad-score model (AUC=0.693, p=0.009) and the clinical model (AUC=0.620, p=0.001) (Supplementary Table 4). On the testing set, the combined model’s sensitivity, specificity, accuracy, precision, recall, F1 score, PPV, and NPV were 0.654, 0.733, 0.704, 0.586, 0.654, 0.618, 0.586, and 0.786, respectively.
Fig. 5A and B display the receiver operating characteristic curves for the Rad-score, clinical, and combined models. Fig. 5C shows the nomogram developed to predict the probability of OPLL progression, integrating Rad-score, preoperative OPLL type, preoperative OPLL volume, and age. The nomogram provides an individualized risk prediction, allowing for practical use in clinical decision-making. Supplementary Fig. 5A and B depict the calibration curves for the combined model on the training and testing sets, respectively. Both curves indicate good calibration, with the bias-corrected line closely following the ideal diagonal line, suggesting a strong agreement between predicted probabilities and observed outcomes. The mean absolute error was 0.012 for the training set and 0.029 for the test set, indicating minimal deviation from ideal calibration. Supplementary Fig. 5C and D present the DCA for the training and testing sets, respectively. The combined model consistently demonstrated a higher net benefit compared to the Rad-score and clinical models across a wide range of threshold probabilities. Fig. 5D and E depict the waterfall plots of predicted probabilities for the training and testing sets of the combined model, respectively. These plots highlight the model’s ability to effectively distinguish between the 2 classes, particularly in the testing set.
4. SHAP Analysis
Fig. 6A–C illustrates the SHAP analysis results to interpret the contribution of each feature to the model’s predictions. Fig. 6A presents the summary plot of SHAP values, indicating the impact of individual features on the model output. The Rad-score had the highest contribution to the model’s prediction, followed by age. Features such as preoperative OPLL type (segmental, mixed, localized, continuous) and preoperative OPLL volume also contributed, but to a lesser extent. Fig. 6B shows the average absolute SHAP values for each feature, providing an overview of their importance. The Rad-score emerged as the most significant predictor, with a mean SHAP value much higher than the other features. Age also showed a considerable impact, while the different OPLL types and volume had relatively lower contributions. Fig. 6C depicts the SHAP dependence plot, illustrating the individual effects of the features on the model output. It shows the interaction effects, with high Rad-score and age values being strong predictors of the outcome, as reflected by positive SHAP values. The plot suggests that both imaging and clinical variables contribute meaningfully to the model’s decision-making process, enhancing the model’s interpretability in predicting OPLL progression.
DISCUSSION
The findings of our study indicate that the combined model, integrating Rad-score and clinical variables, outperformed both the Rad-score-only and clinical models in predicting the postoperative progression of OPLL, achieving an AUC of 0.751 on the test set. This demonstrates the value of combining quantitative imaging biomarkers with clinical information to improve predictive accuracy in OPLL patients. The Rad-score-only model achieved an AUC of 0.693, highlighting the utility of imaging features for risk stratification, while the clinical model alone showed limited predictive power with an AUC of 0.620. The incorporation of clinical variables, such as age, postoperative OPLL type, and preoperative OPLL volume, with the Rad-score significantly enhanced the model’s discriminative ability. These results emphasize the importance of a multimodal approach, integrating Rad-score and clinical data to provide a more comprehensive prediction model, potentially facilitating more informed clinical decision-making and better individualized treatment strategies for OPLL patients.
Our analysis revealed that age, preoperative OPLL type and initial OPLL volume were significant clinical predictors of postoperative progression in OPLL patients. Specifically, our findings indicate that younger patients are more susceptible to OPLL progression, which aligns with previous studies [9,11,20]. This increased risk may be due to higher metabolic activity and the greater mechanical stress exerted on the spine from an increased ROM in younger individuals. We observed that the mixed type of OPLL exhibited the fastest progression, while the localized type showed the slowest progression. This contrasts with the findings of Doi et al. [9], who identified continuous-type OPLL as a risk factor for progression. The discrepancy could be attributed to the differing criteria used to define progression; Doi et al. [9] measured a 2-mm increase in lesion length in any direction, while our study used an annual volume increase of 7.5%. The larger initial volume in continuous-type OPLL might explain its lower progression rate under our criteria. Additionally, our study analyzed the progression patterns of OPLL in each patient, observing that 5 patients progressed from a localized to segmental or mixed type, 42 from segmental to mixed or continuous type, and 30 from mixed to continuous type. This indicates that the classification of OPLL is dynamic and subject to change, which may contribute to the instability of this factor. Our study is the first to find that patients with smaller initial OPLL volumes are less likely to experience OPLL progression, which may be attributed to the fact that smaller OPLL volumes are often associated with the localized type. Surgical techniques may also influence the rate of progression of OPLL. Although several studies have reported that PDF can inhibit OPLL progression by limiting cervical ROM [9,21-23], our results showed no significant difference in progression rates between PDF and LMP. This may be due to selection bias, as patients at higher risk of progression might be more likely to choose the PDF procedure.
Our study utilized 3D OPLL volume analysis, which, compared to 2-dimensional measurements, provides a more comprehensive representation of the entire lesion and allows for precise evaluation of lesion changes. Additionally, for the first time, we applied radiomics methods to extract high-dimensional features from the OPLL lesion, capturing information that is difficult to observe with the naked eye. Our study addresses some of the key challenges in radiomics, such as the high dimensionality of features and the difficulty of interpretation by clinicians. By employing a combination of LASSO and PCA for feature selection and dimensionality reduction, we effectively minimized the risk of overfitting and improved model stability. The use of PCA allowed us to further compress the information into fewer principal components, retaining the most relevant variance while reducing redundancy. Subsequently, the computation of a single Rad-score, derived from these selected principal components, simplified the interpretation of the model’s output, enhancing its clinical applicability. This approach provides a straightforward and interpretable score that can be easily utilized by clinicians for risk assessment, thereby bridging the gap between complex radiomic analysis and practical clinical use. Our findings illustrate the potential of radiomics and machine learning in predicting postoperative outcomes of OPLL surgery, consistent with recent scholarly work in this domain. For instance, Maki et al. [18] developed a model using machine learning and deep learning to predict surgical outcomes in cervical OPLL patients, achieving an accuracy of 71.9% and highlighting the significance of preoperative JOA scores and imaging features. Similarly, Ito et al. [24] demonstrated a model with 74.6% accuracy for overall complications and 91.7% for neurological complications, surpassing traditional LR models. Moreover, Kim et al. [25] found that a machine learning model using an adaptive reinforcement learning algorithm and downsampling outperformed conventional LR in predicting postoperative C5 palsy in cervical OPLL patients, with an AUC of 0.88 compared to 0.69. Although our study achieved an AUC of 0.751 and an accuracy of 70.4%, which are modest compared to some earlier studies, this may be attributed to our focus on a continuous variable outcome, where machine learning models may face challenges with samples near the decision boundary. Our waterfall plots revealed that our model accurately predicts patients with very high or very low risks of OPLL progression, thereby facilitating surgical approach selection and postoperative management. This capability underscores our model’s effectiveness in accurately identifying patients at extreme risks of progression, essential for optimizing surgical and postoperative strategies.
Our study has several limitations that merit attention. Firstly, while it represents the largest sample of OPLL patients with CT follow-up data to date, the single-center, retrospective design of our research poses a significant limitation. Future prospective data collection and external validation at other centers are essential to demonstrate the generalizability of our model. Secondly, our study did not consider a sufficient range of clinical and radiological factors, such as inflammatory markers, the number of OPLL-involved segments, and occupation rate. Including more of these factors in future research could potentially enhance model performance. Thirdly, our radiomics analysis focused solely on the lesion without considering surrounding structures, which may lead to incomplete information extraction.
PDF Links
PubReader
ePub Link
XML Download
Full text via DOI
Full text via PMC
Download Citation
Supplement
Print
