Nomogram Development and Validation for Predicting Postoperative Recurrent Lumbar Disc Herniation Based on Paraspinal Muscle Parameters

Introduction

Intervertebral disc-related low back pain (LBP), comprising 39% of all cases, is mainly attributed to disc herniation (approximately 30%), with other causes being less common.¹ Lumbar disc herniation (LDH) occurs when the fibrous ring of an intervertebral disc ruptures, leading to the protrusion and exposure of the nucleus pulposus in the surrounding fluid.² This leads to inflammation and compression of nerve roots, causing nerve root damage and subsequent symptoms such as lower back and leg pain.³ In recent years, percutaneous endoscopic transforaminal discectomy (PETD) has emerged as a preferred minimally invasive treatment for LDH. PETD offers advantages such as minimal trauma, fewer complications, and quicker recovery.⁴ However, postoperative recurrent lumbar disc herniation (PRLDH) remains a concern, with an incidence ranging from 3–16% according to the literature.^5–7

The spine’s stability depends on the passive subsystem (vertebrae, intervertebral discs, facet joints, spinal ligaments), active subsystem (paraspinal muscles), and neural control subsystem, which are both independent and interconnected.⁸ Therefore, impairment in any of these components can result in reduced spinal function. Crisco et al⁹ conducted load-bearing experiments on cadaveric lumbar vertebrae with muscles removed, retaining only the ligaments. They observed that lumbar stability decreased when subjected to an average load of 88N, whereas the intact lumbar system in a living body could withstand an average load of 2600N. To explore the impact of paraspinal muscles on PRLDH, we reviewed the relevant literature and found that magnetic resonance imaging (MRI) is one of the gold standard methods for assessing whole-body and local skeletal muscle mass, subcutaneous fat tissue and visceral fat tissue, and for LDH patients, MRI can not only assess intervertebral disc degeneration (IVDD), but also observe the relationship between protruding intervertebral disc and nerve root.^10,11 This presents an opportunity to evaluate skeletal muscle mass using MRI. Several studies have demonstrated a correlation between MRI-measured paraspinal muscle fat infiltration, IVDD, and LBP.^12,13 Such infiltration leads to a loss of functional muscle, affecting spinal stability. Therefore, this study aims to investigate the impact of paraspinal functional muscle, assessed via MRI, and other factors on the occurrence of PRLDH after PETD, and to develop and validate the corresponding Nomogram, so that the risk of PRLDH can be more accurately assessed in clinical work. This study strictly adheres to the guidelines of the “Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD)”, with Supplementary Figure S2 providing detailed information.

Materials and Methods

Patients

This retrospective study included 232 patients diagnosed with LDH who were admitted to our hospital between January 2020 and January 2023. The inclusion criteria were as follows: (1) patients met the diagnostic criteria for LDH;⁴ (2) a minimum symptom duration of at least 3 months, with MRI diagnosis consistent with clinical manifestations and signs, and conservative treatment proving ineffective; (3) PETD treatment planned (see Supplementary information S1 for detailed surgical procedure), with significant postoperative symptom relief. Exclusion criteria were as follows: (1) previous lumbar spine surgeries; (2) concurrent spinal conditions such as tumors, tuberculosis, deformities, and fractures affecting spinal structure; (3) long-term use of medications affecting paraspinal muscles after surgery; (4) uncertainty regarding the affected disc; (5) incomplete clinical data; (6) significant cardiovascular or cerebrovascular diseases, or other congenital disorders. PRLDH is characterized by significant alleviation of postoperative neurological symptoms, succeeded by a recurrence of nerve compression symptoms at the same segment and side, confirmed via imaging examinations, around six months post-surgery.¹⁴ The follow-up duration spanned 6 months, primarily conducted through outpatient visits, electronic medical record tracking, and telephone interviews.

Imaging Assessment Method

Paraspinal Muscle

All patients underwent preoperative MRI (3.0T, Siemens, Germany) in a supine position to obtain T2-weighted axial images (see Supplementary information S2 for details). As reported in the literature, the maximum cross-sectional area (CSA) of the paraspinal muscles is located between the L3/4 and L4/5 disc levels, while the largest CSA of the psoas major muscle is at the L4/5 disc level.^15,16 Hence, we selected the paraspinal muscles at the level of the L4/5 intervertebral disc. The T2-weighted images were imported into Image J software. Two experienced radiologists (reader 1, with 8 years of MRI diagnostic experience; reader 2, with 4 years of MRI diagnostic experience) independently measured the preoperative contralateral and ipsilateral multifidus CSA, erector spinae CSA, and psoas major CSA. The fat CSA of the contralateral and ipsilateral multifidus, erector spinae, and psoas major were measured. Mean values were used for statistical analysis (Figure 1). The functional cross-sectional area (FCSA) of the paraspinal muscles was calculated as FCSA = Paraspinal Muscle CSA - Paraspinal Muscle Fat CSA. Additionally, to control for inter-individual size differences, the FCSA was transformed into the standardized functional cross-sectional area (SFCSA), given by SFCSA = Paraspinal Muscle FCSA / Corresponding Vertebral Body CSA × 100%. To further normalize for different heights among patients, the paraspinal muscles’ skeletal muscle index (SMI) was calculated as SMI = Paraspinal Muscle FCSA (mm²) ÷ Height2 (m²).¹⁷

Figure 1 L4/5 paraspinal muscles at the lumbar interspace. (A) T2-weighted images (MF: Multifidus muscle; ES: Erector spinae; PS: Psoas major). (B) Computation of paraspinal muscle cross-sectional area and fatty infiltration area (a: Multifidus muscle; b: Erector spinae; c: Psoas major).

Vertebral Bodies and Facet Joints

Using T2-weighted sagittal MRI images, we measured the IA between the upper and lower vertebral bodies of the affected disc, and by subtracting the IAs of the upper and lower vertebral bodies, we obtained the intervertebral disc angle (IDA) (Figure 2). Additionally, using CT axial images (MX8000, MACONI, USA), we measured the facet joint angle (FJA) of the affected disc. The FJA was determined as the angle between the line connecting the anterior and posterior points of the facet joint and the tangent to the posterior margin of the vertebral body at the lower level of the intervertebral disc (Figure 2).

Figure 2 (A) T2-weighted MRI. γ: Upper vertebral endplate inclination angle (IA); Ω: Lower vertebral endplate IA; Intervertebral disc angle (IDA) = γ - Ω. (B) c: Tangent to the posterior edge of the vertebral body; ab: Line connecting the facet joints (FJ) anterior-medial point and the FJ’s posterior-lateral point; Right facet joint angle (FJA): Angle α formed by points a and c; Left FJA: Angle β formed by points b and c.

Risk Factors

The included risk factors were as follows: General characteristics: age, gender, height, body mass index (BMI), diabetes, hypertension, smoking, drinking, preoperative visual analog scale (VAS) score, postoperative VAS score, educational level, occupation, course of disease, hospital stay, and duration of surgery. Imaging features: Pfirrmann grade, Modic changes, IDA, left and right FJA, contralesional multifidus muscle’s SFCSA and SMI, ipsilesional multifidus muscle’s SFCSA and SMI, contralesional erector spinae muscle’s SFCSA and SMI, ipsilesional erector spinae muscle’s SFCSA and SMI, contralesional psoas major muscle’s SFCSA and SMI, ipsilesional psoas major muscle’s SFCSA and SMI, herniation type, and herniation segments.

Statistical Analysis

Statistical analyses were conducted using R Studio software (v4.2.3, http://www.rproject.org/). The “irr” package was employed to compute the intraclass correlation coefficient (ICC), with results ≥0.75 indicating good consistency. The “glmnet” package performed Least Absolute Shrinkage and Selection Operator (LASSO) regression for the selection of predictive factors for PRLDH occurrence. Multiple logistic regression was carried out using the “rms” package. The “nomogramFormula” package was utilized to create the Nomogram. Calibration curves and the Hosmer-Lemeshow goodness-of-fit test were generated using the “rms” package. Decision Curve Analysis (DCA) and Clinical Impact Curve (CIC) were plotted using the “rmda” package. Statistical significance was set at P < 0.05. All modifications to the typeface and font size of the figures were made using Adobe Illustrator. Additional details are available in Supplementary Information S3.¹⁸

Results

Comparison of Clinical Data Between Two Groups

This study included 232 patients, comprising 140 males and 92 females, with an age range of 21 to 83 years (mean age: 52.94±13.88 years). Among these patients, 26 experienced PRLDH, resulting in a recurrence rate of 11.21%. Clinical information for recurrent and non-recurrent patients is compared in Supplementary Tables S1 and S2. Excellent consistency in data measurements between the two observers was observed (0.847–0.922, Supplementary Table S3). Using the “caret” package in R, random sampling was performed with the createDataPartition function, dividing collected cases into training (n=140) and validation (n=92) sets in a 6:4 ratio (Table 1).

Table 1 Characteristics of Training and Validation Groups

Model Development

To prevent overfitting, each feature requires participation from at least 10–15 patients in model construction.^19,20 Given the training group’s 140 patients, the maximum feature limit is set at 14. Utilized Lasso regression for variable selection. The model exhibited optimal performance (lambda = 0.056) with four predictors (Figure 3). Diabetes, Modic changes, ipsilesional multifidus SMI, and ipsilesional erector spinae SMI were incorporated into a multivariable logistic regression model. Ultimately, diabetes, Modic changes, and ipsilesional multifidus SMI (P < 0.05) emerged as independent predictors for PRLDH post-PETD (Table 2). Developed a Nomogram model with these predictors (Figure 3) and made it available at https://sofarnomogram.shinyapps.io/DynNo_PRLDH/.

Table 2 Multivariate Logistic Regression Analysis

Figure 3 (A) Characterizing the Variations of LASSO Regression Coefficients. (B) LASSO regression selects the optimal parameter lambda through cross-validation. The dashed line on the right represents lambda values with average error within ±1 standard deviation, indicating improved model performance. (C) Developed Nomogram based on multifactorial logistic regression analysis.

Model Validation

Through 1000 Bootstrap resampling iterations, the training group’s AUC was 0.921 (95% CI 0.872–0.970), and the validation group’s AUC was 0.900 (95% CI 0.828–0.972) (Figure 4). Demonstrating good predictive performance, the model is detailed in Table 3 and Supplementary Figure S1. Calibration curves for both groups exhibited well-fitted relationships between observed and predicted values (Figure 4). The Hosmer-Lemeshow goodness-of-fit test showed good fit in the training group with χ² = 5.638, P= 0.688, and in the validation group with χ² = 6.259, P= 0.618.

Table 3 Diagnostic Performance of the Nomogram and Ipsilesional Multifidus SMI Model for Predicting PRLDH

Figure 4 (A and B) The ROC curves of the training and validation groups exhibit AUC values of 0.921 (95% CI 0.872–0.970) and 0.900 (95% CI 0.828–0.972), respectively. (C and D) Training and validation set calibration curves.

Clinical Utility of the Predictive Model

According to DCA curve results in the training and validation groups, the predictive model yields the maximum clinical net benefit within risk thresholds of 0.02–1.00 and 0.02–0.80. It provides significant additional clinical net benefit in predicting PRLDH occurrence (Figure 5). Risk stratification for 1000 cases using CIC indicates (Figure 5) that within the threshold probability range, the predicted number of patients experiencing PRLDH consistently exceeds the actual occurrences.

Figure 5 (A and B) The prediction probability and the number of patients in the training set and the validation set of the nomogram, the dashed line indicates the cutoff value when the maximum Youden index is reached. (C and D) Decision Curve Analysis (DCA) curves for the training and validation sets. (E and F) Clinical Impact Curve (CIC) for the training and validation sets.

Discussion

In our study, we identified diabetes, Modic changes, and ipsilesional multifidus SMI as independent predictors for PRLDH. The model demonstrated good predictive performance in both the training and validation sets. Meanwhile, the model’s performance was unaffected by age (≥60 years) (AUC: 0.914 vs 0.950; P = 0.491) or gender (AUC: 0.915 vs 0.972; P = 0.146) (Supplementary Table S4). The observed incidence of PRLDH in our study was 11.21%, aligning with the reported range in the literature (3–16%). In the current study, only a few have utilized paraspinal muscle parameters to predict PRLDH, primarily focusing on fat infiltration grades.²¹ Quantitative studies on paraspinal muscles are as follows: Choi et al²² conducted a retrospective analysis of 49 patients who underwent discectomy and partial laminectomy, revealing larger psoas major muscle parameters in PRLDH patients compared to those No-PRLDH. However, a larger study by Kong et al⁶ (654 patients undergoing PETD) found smaller SFCSA at the L4/5 level in PRLDH patients (excluding the psoas major muscle), consistent with our findings. Nevertheless, previous studies have not explored the impact of different paraspinal muscles or asymmetry in paraspinal muscles on PRLDH occurrence. Our study demonstrates that paraspinal muscle SMI parameters can predict PRLDH occurrence, and we found no significant differences in contralesional and ipsilesional paraspinal muscle parameters among recurrent patients. However, most paraspinal muscle parameters in recurrent patients were significantly lower than those in non-recurrent patients, indicating pronounced paraspinal muscle degeneration (fat infiltration, fibrosis, and atrophy) in recurrent cases. Ipsilesional multifidus SMI may act as a protective factor against PRLDH. Its mechanism may involve nerve root injury resulting from disc herniation, leading to structural changes in ipsilateral paraspinal muscles, such as increased fat infiltration and type I to II fibers conversion. These changes may contribute to lower back pain and declining spinal function.²³ Therefore, prompt relief of nerve root compression is crucial for preserving paraspinal muscle mass, which not only alleviates patient discomfort but also improves prognosis.

Panjabi’s⁸ groundbreaking elucidation of the “three-system” theory and the concept of the “neutral zone” in spinal biomechanics has significantly advanced our understanding of the complexity and functional integration of the spinal structure. This has propelled research into other components of the spine, such as vertebral bodies, ligaments, and facet joints, enabling a better comprehension of the intricacies of lumbar spine disorders. Among these structures, paraspinal muscles play a crucial role in maintaining spinal stability, particularly in the “neutral zone” and lumbar stability. As a local stabilizer, the multifidus muscle accounts for at least 2/3 of maintaining the stability of the “neutral zone”, and when the muscle degenerates, it will reduce the ability to control the neutral zone, and pain improvement does not necessarily mean multifidus muscle function recovery, but multifidus muscle atrophy may be related to recurrent LBP.^24,25 According to Zhu et al’s²⁶ study, PRLDH is more likely to occur when the multifidus muscle at the herniation site exhibits moderate to severe fatty infiltration. Studies have shown that the physiological cross-sectional area (PCSA) of the lumbar multifidus muscle is larger than that of the erector spinae, enabling it to generate more force.^27,28 This is attributed to the muscle fiber composition (I, IIA, and IIX types) of the paraspinal muscles. Agten et al’s²⁹ study indicates that, compared to the erector spinae, the multifidus has significantly larger type I muscle fibers, highlighting its greater role in maintaining spinal stability. Meanwhile, multifidus muscles are crucial for stabilizing and controlling spinal movement.²⁵

This study quantified paraspinal muscle FCSA and further calculated paraspinal muscle SFCSA and SMI to account for inter-individual differences. While muscle volume typically correlates with patient height, degeneration of the vertebral body can lead to a decrease in the CSA/vertebral body CSA ratio due to diffuse proliferation, potentially compromising the reliability of SFCSA measurements.³⁰

In our study, Modic changes and diabetes were also considered as predictive factors. The vertebral endplate consists of a bony endplate and a cartilaginous endplate, where the bony endplate is formed by ossification of the vertebral surface epiphyseal plate, while the cartilaginous endplate consists of a central thin layer of translucent cartilage of the vertebral endplate. Modic changes are defined based on abnormal MRI alterations of the vertebral endplate and subchondral bone. Previous research has indicated that lumbar Modic changes may contribute to lower back pain.³¹ Modic changes may imply a state of intense inflammation. According to Luo et al’s³² research, the use of exosomes derived from normal cartilage endplate stem cells (N-Exos) and degenerated cartilage endplate stem cells (D-Exos) in rats demonstrated their impact on apoptosis in nucleus pulposus cells and IVDD. They found that N-Exos were more effective in activating autophagy to suppress IVDD compared to D-Exos. Additionally, normal cartilage endplate stem cells inhibited apoptosis in nucleus pulposus cells by suppressing the inflammatory response. Both cartilage endplate degeneration and nucleus pulposus cell apoptosis are significant triggers of IVDD. When herniated material is exposed to body fluids, it triggers an inflammatory response, leading to the deterioration of the intervertebral disc microenvironment. Moreover, surgical procedures can damage the vertebral endplate, further exacerbating the microenvironment deterioration, thereby increasing the risk of PRLDH occurrence.

A meta-analysis by Luo et al³³ involving 58 cohort studies revealed a significant association between PRLDH and diabetes. This correlation may be linked to diabetes’ impact on the body’s connective tissues. Kakadiya et al³⁴ conducted a study on intervertebral disc tissues of diabetic patients using histopathological and immunofluorescence analysis, and observed reduced proteoglycan synthesis, increased aggrecan cleavage, and elevated expression of ADAMTS4/5 proteins (involved in proteoglycan breakdown within the intervertebral disc). These findings indicate that hyperglycemia accelerates apoptosis in intervertebral disc tissues. Meanwhile, an animal study demonstrated that high glucose increases stress-induced senescence (p16-pRB) pathway-related protein expression in a dose- and time-dependent manner, leading to senescence of nucleus pulposus cells in adult rats.³⁵ In conclusion, our developed Nomogram predictive model indicates that a smaller SMI of ipsilesional erector spinae or multifidus muscles, higher Modic change grade, and the presence of diabetes significantly increase the risk of PRLDH occurrence.

This study has several limitations. Firstly, this study is a single-center retrospective study. Although we performed internal validation using the training group, validation group, and enhanced Bootstrap method, additional clinical cohorts are needed for external validation. Secondly, we only included three predictive variables, while many other factors such as physical activity and inflammatory responses were not considered. Future research should address these limitations and further explore the impact of additional relevant factors.

Conclusion

This study found differences in paraspinal muscle parameters between recurrent and non-recurrent patients, with no significant asymmetry observed in the contralesional and ipsilesional parameters of recurrent patients. Moreover, ipsilesional multifidus SMI was identified as a protective factor for PRLDH occurrence. Finally, the Nomogram predictive model, developed based on paraspinal muscle parameters in this study, demonstrates excellent predictive capability, facilitating risk assessment for the occurrence of PRLDH after PETD.