Between January 2015 and January 2022, a total of 876 patients with single-segment OVCFs who underwent unilateral PKP in our center were enrolled. Patients in the study satisfied the following inclusion criteria: 1) aged > 55 years, 2) first-time PKP procedure, 3) single segment vertebral fracture of T4-L5, 4) thoracolumbar MRI suggesting bone marrow edema (hypointense T1 and hyperintense T2), 5) diagnosed with osteoporosis by dual energy X-ray absorptiometry (T score ≤ -2.5), and 6) obvious thoracolumbar and back pain with a score of six or more on the visual analogue scale (VAS) and limited physical activity. The following were the exclusion criteria: 1) patients with other medical conditions that might cause thoracolumbar and back pain, such as infections, malignant tumors, and a history of spinal surgery, 2) patients whose spinal canal invaded, 3) unable to tolerate surgery, such as when severe cardiopulmonary comorbidity is present, and 4) patients with insufficient follow-up data.

In the end, 731 strictly screened OVCF patients were recruited. Patients between January 2017 and December 2020 were involved in a training cohort (n = 548) and those between January 2021 and January 2022 were included in a validation cohort (n = 183). Figure 1 shows the flowchart illustrating patient selection and model design.

The baseline data included age, gender, body mass index (BMI), comorbidities (hypertension, diabetes), smoking status, BMD, fracture position, as well as preoperative VAS score and oswestry disability index (ODI). Radiological parameters, including vertebral height loss, Cobb angle (°) (the angle between the superior and inferior endplate of the fractured level), presence of IVC, and TLF injury were obtained from preoperative X-ray, CT, and MRI imaging.

Patients were placed in the prone position for PKP while being given local or general anaesthetic. A bone puncture trocar was placed at the broken level via the lateral pedicle edge and gradually advanced through the pedicle into the vertebral body under C-arm fluoroscopic guidance. The vertebral body was then carefully injected with polymethylmethacrylate using an inflated bone tamp. If cement spilled into extraosseous structures or veins, reached the cortical margin of the vertebral body, or both, the injection was terminated. On the first day following PKP, ambulation was suggested for all patients. Alendronate (70 mg/week), calcitriol (0.25–0.5 g/day), and calcium tablets (600–1200 mg/day) were given postoperatively. Non-steroidal anti-inflammatory medicines were administered to patients with inadequate pain relief or obvious pain that interfered with daily life.

Patients were followed up and the VAS scores were recorded at four different time intervals, including 1, 3, 7 and 30 days after surgery. RBP was defined as the VAS score ≥ 4 at both 3 and 30 days based on previous experience and literature review due to the lack of existing guidelines (9, 17).

CT scans were acquired from a 64-slice spiral CT scanner (Aquilion Prime Model TSX-303A, Toshiba Medical Systems Corp., Tokyo, Japan) using the following scanning parameters: tube voltage 120-130 kV, variable mAs to achieve a noise index of 25, rotation time 0.5 or 0.75 s, and matrix 512 × 512. In all cases, axial views with a slice thickness of 3 mm were assessed initially, followed by sagittal and coronal reformations using the inhouse CT postprocessing software.

The radiomics workflow is shown in Figure 2 , which mainly includes image segmentation, feature extraction, and feature selection. Vertebral segmentation was performed separately by two radiologists with more than 10 years of clinical experience using 3D Slicer software (v4.10.2) to manually delineate an ROI containing the entire vertebral body slice by slice on axial images. The radiomics features were extracted with the open-source Pyradiomics package. There were 6 image types (original, exponential, gradient, logarithm, square, and squareroot) and 7 feature classes including shape, first-order statistics (first-order), gray-level co-occurrence matrix (glcm), gray-level dependence matrix (gldm), gray-level run length matrix (glrlm), gray-level size zone matrix (glszm), and neighbouring grey tone difference matrix (ngtdm) adopted for each sequence (18). We employed wavelet decomposition filtering after image reconstruction.

All features extracted from each segmented image were normalized using the standard score (z-score) transformation. To avoid overfitting and improve feature repeatability, three steps were taken in the feature selection process. First, the intraclass correlation coefficient (ICC), with an ICC > 0.9 considered reliable, was used to evaluate discrepancies between the features delineated by two radiologists. Each stable feature was then compared between patient cohorts with and without RBP using a Student t-test. Finally, the features with false discovery rate (FDR)-adjusted P values < 0.05 using Benjamini–Hochberg correction were then subjected to least absolute shrinkage and selection operator (LASSO) regression analysis. Each rad-score is constructed by weighting the coefficients of selected radiomics features based on the output of the LASSO regression.

The rad-scores and those RBP predictors initially screened by univariate analyses were included in the subsequent multivariate analysis to determine the independent correlates of RBP. These independent variables were then used to create a nomogram for preoperatively identifying whether an OVCF patient will suffer RBP. Prior to being externally evaluated in the validation cohort, its performance was first assessed in the training cohort.

In radiomics, feature extraction, normalization, selection, and rad-score creation were run under Python 3.7.1. Univariate and multivariate Logistic regression analyses were applied to screen for the independent correlates of RBP. The rad-scores plus those RBP predictors that had been first assessed by univariate studies were involved in the subsequent multivariate analysis to determine the odds ratio (OR) and corresponding 95% confidence interval (CI) of each predictor. The nomogram was established according to these correlates and then validated in the training and validation cohorts. Bootstraps analyses were utilized to evaluate the unbiased performance of the model. We evaluated the discrimination and calibration of the nomogram by the Harrell’s concordance index (C-index) combined with the receiver operating characteristic (ROC) curve, and the calibration curve along with the Hosmer-Lemeshow (HL) test. The clinical utility was assessed using a decision curve analysis (DCA), which quantified the net benefit of each decision strategy at each threshold probability. SPSS (version 20.0), MedCalc software (version 19.2.1), and the R package (version 4.2.1) were used for all statistical analyses not related to radiomics.

Finally, 731 OVCF patients were enrolled (548 patients in the training cohort for model development and 183 in the validation cohort for evaluating model performance). RBP accounted for 11.2% (81/731) in the overall cohort while 10.8% (59/548) and 12.0% (22/183) in the training and validation cohorts, respectively. There were no statistically significant differences in baseline data and preoperative radiological parameters between the training and validation cohorts (all P values > 0.05). The details of two the cohorts are listed in Supplemental Table 1 .

In the training cohort, 1316 radiomics features for each segmented image were extracted and normalized. Of these, 1077 (81.8%) features with favorable reproducibility (ICC > 0.9) were selected for subsequent t-test screening, resulting in 85 initially screened radiomics features. By constructing a penalty function, 8 radiomics features were chosen by determining the best penalty regularization parameter (λ) with the 1-standard error for the minimum criteria in the LASSO model ( Figure 3 ). The rad-score of each patient was established based on the weights of these 8 features and their corresponding coefficients. The calculation of rad-score is shown in Supplemental Appendix 1 .

The results of univariate and multivariate analyses to determine the independent correlates of RBP are summarized in Table 1 . Rad-score plus BMD, IVC, and TLF injury were incorporated into the multivariate analysis after being filtered by the univariate analyses. The result revealed that RBP was independently associated with BMD, IVC, TLF injury, and rad-score (all P < 0.05). RBP was more prevalent in OVCF patients who had IVC, TLF injury, higher rad-score, and lower BMD (OR: 3.220, 11.377, 1.062, and 0.082, respectively).

Using these selected variables, a radiomics nomogram was developed to preoperatively visualize the likelihood of RBP in OVCF patients. As indicated in Figure 4 , RBP risk was assessed by accumulating the points corresponding to each variable. For instance, the BMD and rad-score of a 70-year-old female OVCF patient with TLF injury are -3.1 and 0.9, respectively. The corresponding scores are as follows: 0 points for IVC, 60 points for TLF injury, 30 points for BMD, and 60 points for rad-score. She has about a 60% chance of encountering RBP based on her overall score of 150.

The discrimination, calibration, and clinical utility for the developed nomogram were first evaluated in the training cohort. Following 1,000 bootstrapping, the stability of the model was tested to adjust the overfitting deviation. An adjusted C-index of 0.936 and a nonsignificant P value of the HL test (χ² = 5.796, P = 0.670) indicated good discrimination and calibration power, as indicated in Figures 5A, B . Moreover, the result of DCA, which is displayed in Figure 5C , demonstrated that stratifying the RBP of OVCF patients using the nomogram would obtain clinical net benefits.

Subsequently, this model was independently tested in the validation cohort and found that it still revealed satisfied discrimination (C-index = 0.787) ( Figure 5D ) and calibration (χ² = 8.327, P = 0.238) ( Figure 5E ). The DCA plot in Figure 5F suggested that although the clinical net benefit of the application in the validation cohort to predict the RBP risk was lower than in the training cohort, it still had good potential for clinical utility.

To further highlight the complementary role of radiomics in prognostic prediction, we compared the predictive performance of the rad-score itself, and the nomograms with and without the rad-score. As shown in Figure 6 , the prediction accuracies of the rad-score itself and the nomogram that only incorporated regular features were similar (AUC: 0.797 and 0.810, respectively), and both were significantly lower than the radiomics nomogram (AUC: 0.936, both P < 0.05). It implied that the model combining human-recognizable and computer-extracted features was more accurate for prognosis prediction.