Obesity Induced Intradiscal Stress Amplification Following Lumbar Foraminoplasty: Implications for Reherniation Risk

2.1. Geometric Parameters and Material Properties

The lumbar spine model was developed from computed tomography (CT) images of a 55-year-old Indonesian man. He had a body weight of 66 kg and a height of 168 cm. His bone condition was considered normal based on both clinical evaluation and radiological examination.

The CT scan data featured high-resolution images of the lumbar spine with a slice thickness of 1.0 mm. A single CT-derived anatomy was used because the objective of this study was to perform a controlled within-model comparison of BMI-derived axial loading rather than to estimate population-level variability. The same geometry, material properties, contact definitions, surgical configuration, and boundary conditions were maintained across all five loading scenarios. Consequently, differences in the predicted stress responses could be evaluated primarily as a function of axial load magnitude. Similar single-subject CT-based finite element approaches have been used for controlled comparative analyses of lumbar biomechanics [18, 20]. Data collection was performed with the volunteer’s informed consent. The subject was recruited through the Department of Spinal Surgery at Dr. Kariadi General Hospital, Semarang, Indonesia. The study protocol was approved by the Ethical Committee of Medical Faculty of Medicine, Diponegoro University, Semarang, Indonesia (Ethical Approval No. 108/EC/KEPK/FK-UNDIP/IV/2022). The study was conducted in accordance with the principles of the Declaration of Helsinki.

The imaging data were then processed using Mimics 21.0 (Materialise, Belgium) to segment the L4–L5 vertebrae and generate the initial three-dimensional geometry as the base model. The STL file derived from this initial geometry was subsequently refined and prepared for analysis in ANSYS Workbench 2023 R2 (Ansys, USA) prior to testing.

To distinguish between cortical and cancellous bone, each vertebral body was modeled as a cancellous core surrounded by a cortical shell. The cortical shell was generated along the outer vertebral surface with a variable thickness of 0.5–1.0 mm, while the remaining internal volume was assigned as cancellous bone. This modeling approach allowed the vertebrae to preserve the anatomical distinction between the stiffer outer cortical layer and the less dense internal cancellous region.

Maintaining consistency between the material assignment in a finite element model and the actual properties of spinal tissues is considered important for achieving a more reliable biomechanical representation. Accordingly, the material properties of the modeled spinal structures were assigned based on previously reported finite element parameters [21]. The material properties used for each spinal structure, including elastic modulus and Poisson’s ratio, are presented in Table 1 to support model fidelity and to provide a consistent reference for future simulations.

To make the model biomechanically closer to the real spinal condition, seven ligament structures were also included: the anterior longitudinal ligament (ALL), posterior longitudinal ligament (PLL), interspinous ligament (ISL), supraspinous ligament (SSL), ligamentum flavum (LF), intertransverse ligament (ITL), and capsular ligament (CL). These ligaments were modeled as nonlinear spring elements with tension-only behavior, which follows their physiological role. Most of the ligament material properties were adopted from a previously published finite element model [21]. A complete summary of the modeling parameters is provided in Table 2.

2.2. Model PTED Foraminoplasty

The finite element model was developed and focused on a single functional spinal unit (FSU), the L4–L5 segment. This segment was chosen because it is one of the lumbar levels most often involved in disc herniation and stenosis, making it highly relevant from a clinical perspective. The postoperative geometry of the L4–L5 FSU was reconstructed from anatomical references. The goal was to create a post-foraminoplasty model that could approximate the actual postoperative condition more closely. A defined resection was also made at the apex of the L5 superior articular process (SAP). This was intended to represent extensive decompression, a procedure often required when the neural foramen is severely narrowed.

The resection configuration used in this study was determined with reference to a previously published finite element framework [16] which categorized foraminoplasty resections into five groups (A–E). Group A served as the intact model. Group B involved ventral resection of the superior articular process without exposing the articular surface. Group C, the apex resection group, included excision of the SAP apex with partial exposure of the articular surface and removal of the adjacent capsular ligament. Groups D and E involved progressively more extensive resections, incorporating removal of the isthmus and lateral recess, and in Group D, even part of the L5 pedicle notch.

Based on previous findings, Group B resections produced negligible changes in intervertebral disc stress and range of motion (ROM), rendering them unsuitable for evaluating stress redistribution. Conversely, Groups D and E induced substantial increases in intervertebral disc stress (IDS) and kinematic instability, making them less clinically representative and potentially unsuitable for stable finite element simulations. Group C was therefore selected for the present study, as it represents a biomechanical threshold in which the facet joint is partially compromised, allowing for meaningful stress redistribution without introducing excessive instability. This configuration also reflects a clinically realistic decompression extent frequently performed during moderate foraminoplasty. The geometric characteristics of the adopted resection, including articular surface exposure, are illustrated in Fig. (1).

2.3. Boundary Condition

To replicate physiological loading conditions of the lumbar spine, the finite element model of the L4–L5 segment was subjected to a compressive force of 500 N to simulate the upper body load, combined with an applied moment of 10 N·m to induce representative spinal motions, including flexion, extension, lateral bending, and axial rotation [21]. The compressive load was applied to the superior surface of the L4 vertebra, while the inferior surface of the L5 vertebra was fixed in all degrees of freedom to provide a stable boundary condition. The applied moment was imposed on the superior region of L4 according to the corresponding direction of each simulated motion. These loading scenarios reflect typical biomechanical demands experienced during routine activities of daily living. The resulting segmental ROM and von Mises stress distribution within the intervertebral disc were computed to evaluate the mechanical behavior of the spinal segment.

To enhance physiological relevance, the boundary conditions were further adapted using BMI-derived axial compressive loading scenarios. The compressive loads were scaled using factors derived from the estimated body mass of each BMI category relative to the normal-weight condition. Accordingly, the applied axial loads were 394.65 N, 500 N, 615.89 N, 665.93 N, and 794.98 N, representing underweight, normal weight, at-risk of obesity, obesity Class I, and obesity Class II conditions, respectively [22]. The boundary conditions and anatomical configuration of the lumbar segment are illustrated in Fig. (2).

2.4. Mesh Study

The solid components of the L4–L5 finite element model were discretized using ten-node quadratic tetrahedral elements. Mesh convergence was evaluated by progressively refining the element size and monitoring both maximum von Mises stress in the annulus fibrosus and total deformation. Maximum annulus fibrosus stress was selected as the principal convergence parameter because intradiscal stress was the primary outcome of the present study, whereas total deformation was used to evaluate the stability of the global structural response.

The predicted maximum annulus fibrosus stress varied at coarser element sizes but stabilized after the element size reached 1.0 mm. The maximum stress values obtained at element sizes of 1.0, 0.9, and 0.8 mm were 1.812, 1.800, and 1.810 MPa, respectively. The successive differences were approximately 0.66% between 1.0 and 0.9 mm and 0.56% between 0.9 and 0.8 mm. These differences were below the predefined convergence criterion of 1%, indicating that further mesh refinement produced only negligible changes in the predicted annulus fibrosus stress.

Total deformation was also assessed as a secondary convergence parameter. At an element size of 1.0 mm, the complete L4–L5 model contained 588,628 elements. Refinement to 0.9 mm increased the total number of elements to 800,803, representing an increase of approximately 36%, while the resulting total deformation differed by only 0.2%. Therefore, the 1.0-mm mesh was adopted for the subsequent simulations because both the local annulus fibrosus stress and global deformation responses had stabilized, while further refinement substantially increased the computational cost. The convergence behavior of the finite element model across different element sizes is illustrated in Fig. (3).

2.5. Validation Method

The finite element model was validated by comparing the predicted range of motion (ROM) of the L4–L5 segment with previously reported in vitro experimental data [23]. The simulations were performed under flexion–extension, lateral bending, and axial rotation loading conditions. To improve the physiological relevance of the model, the boundary conditions and external loading configuration were adopted from a previously published finite element method [21]. These loading and constraint conditions were used for model implementation and were not treated as a separate quantitative validation dataset.

Quantitative agreement was assessed using the absolute deviation and percentage deviation between the predicted ROM and the corresponding experimental mean. The absolute deviations were also expressed relative to the reported experimental standard deviations to account for inter-specimen variability. In addition, each predicted value was evaluated against the corresponding experimental mean ± standard deviation interval. The mean absolute error (MAE) across the three loading directions was calculated as an overall measure of prediction error. Because the present study employed a deterministic, single-geometry finite element model with one numerical prediction for each motion condition, the statistical assessment was descriptive rather than inferential; therefore, hypothesis testing and p-values were not applied.

The mean absolute error (MAE) was calculated to summarize the overall magnitude of the differences between the predicted and experimental mean ROM values across the three motion conditions:

where P_i represents the predicted ROM, R_i represents the corresponding experimental mean ROM, and nis the number of motion conditions evaluated (n = 3).

2.6. Quantitative Comparison of Stress Responses

Maximum von Mises stress was compared descriptively across the five BMI-derived axial loading scenarios for each motion direction and disc component. The 394.65 N condition was used as the reference. Absolute stress differences were calculated by subtracting the reference stress from the stress under each higher loading condition, while percentage changes were calculated relative to the reference value. Because the analysis was based on a deterministic single-geometry finite element model with one numerical output for each condition, inferential statistical tests and p-values were not applied.

3.1. Model Validation

The quantitative validation results are presented in Fig. (4) and Table 4. The predicted ROM values were 14.70° for flexion–extension, 10.20° for lateral bending, and 2.50° for axial rotation. Compared with the reference experimental mean values [23], the corresponding absolute deviations were 0.10°, 2.00°, and 1.20°, representing percentage deviations of 0.68%, 16.39%, and 32.43%, respectively. The overall MAE across the three motion conditions was 1.10°.

Despite the differences relative to the experimental means, all predicted ROM values remained within the corresponding experimental mean ± standard deviation intervals. The predicted flexion–extension value was within the experimental interval of 12.70–16.90°, lateral bending was within 9.95–14.45°, and axial rotation was within 2.20–5.20°. When the absolute differences were expressed relative to the corresponding experimental standard deviations, the deviations were 0.05, 0.89, and 0.80 standard deviations for flexion–extension, lateral bending, and axial rotation, respectively.

Axial rotation exhibited the highest percentage deviation. However, this result should be interpreted in conjunction with the absolute deviation and experimental variability. Because the experimental mean ROM for axial rotation was relatively small, an absolute difference of 1.20° produced a comparatively large percentage deviation. Nevertheless, the predicted axial-rotation value remained within one experimental standard deviation of the reference mean. The remaining discrepancies may be associated with subject-specific anatomical characteristics, differences in vertebral and facet geometry, the use of simplified isotropic linear-elastic material properties, the representation of ligaments using spring elements, variations in contact definitions, differences in boundary and loading conditions, and the inherent inter-specimen variability of experimental spinal measurements. Overall, the agreement with the experimentally reported ROM ranges supports the use of the model for the subsequent comparative biomechanical analyses.

3.2. Von Mises Stress on Vertebral Disk L4-L5

As shown in Fig. (5), maximum von Mises stress in the nucleus pulposus generally increased with greater BMI-derived axial compressive loading. Flexion consistently produced the highest nucleus pulposus stress across all five loading conditions, increasing from 0.41650 MPa at 394.65 N to 0.58430 MPa at 794.98 N. This corresponded to an absolute increase of 0.16780 MPa and a relative increase of 40.29%.

From the 394.65-N to the 794.98-N condition, nucleus pulposus stress increased by 246.61% during extension, 80.90% during left lateral bending, 90.08% during right lateral bending, 42.37% during left axial rotation, and 81.81% during right axial rotation. The large percentage increase during extension should be interpreted together with its low reference stress of 0.07771 MPa; the corresponding absolute increase was 0.19164 MPa, and the final extension stress remained lower than the flexion stress.

The response was not strictly monotonic at every adjacent loading step. A small decrease occurred during left axial rotation between 394.65 and 500 N and during right lateral bending between 615.89 and 665.93 N. Therefore, the results indicate an overall load-dependent increase rather than a uniformly progressive increase under every condition.

As shown in Fig. (6), maximum von Mises stress in the annulus fibrosus also demonstrated an overall increase with greater BMI-derived axial loading, although the response was not strictly monotonic and the motion producing the maximum stress varied among loading levels. During flexion, annulus fibrosus stress increased from 1.94500 MPa at 394.65 N to 2.44200 MPa at 794.98 N, corresponding to an absolute increase of 0.49700 MPa and a relative increase of 25.55%.

From the 394.65 N to the 794.98 N condition, annulus fibrosus stress increased by 21.99% during extension, 41.28% during left lateral bending, 17.26% during right lateral bending, 33.03% during left axial rotation, and 49.61% during right axial rotation. The overall maximum annulus fibrosus stress was 2.66130 MPa during left lateral bending at 794.98 N.

The dominant motion varied according to loading level. Right lateral bending produced the highest annulus fibrosus stress at 394.65 N, flexion at 500 and 615.89 N, right lateral bending at 665.93 N, and left lateral bending at 794.98 N. These findings indicate that annulus fibrosus stress was influenced by both axial load magnitude and motion direction.

Across the evaluated conditions, the magnitude and motion producing peak stress differed between the nucleus pulposus and annulus fibrosus.

3.3. Regional Stress Distribution

The stress contours in Fig. (7) showed that nucleus pulposus stress during flexion was predominantly concentrated in the posterior region. Across the evaluated loading conditions, the principal stress concentration remained generally posterior, although its magnitude varied with axial load.

Figure 8 showed two principal stress-concentration regions in the annulus fibrosus, located in the anterolateral and posterolateral regions. The relative prominence of these regions varied according to motion direction and axial loading condition.

4.1. Influence of Axial Load Magnitude and Motion Direction

Overall, greater BMI-derived axial loading increased mechanical demand on both disc components, but the nucleus pulposus and annulus fibrosus exhibited different motion-dependent responses. The hypothesis was therefore partially supported: stress generally increased in both tissues and flexion consistently produced the highest nucleus pulposus stress, whereas the maximum annulus fibrosus stress varied among motion directions.

The simulations showed that increasing BMI-derived axial loading generally increased intradiscal stress, although the magnitude of the response differed among motion directions and was not uniform at every loading increment. Between 394.65 and 794.98 N, nucleus pulposus stress increased by 40.29–246.61%, whereas annulus fibrosus stress increased by 17.26–49.61%. Several small decreases occurred between adjacent loading conditions, but these variations did not alter the overall pattern of greater stress under the highest axial load.

The largest relative increase in nucleus pulposus stress occurred during extension, rising from 0.07771 to 0.26935 MPa. This percentage was amplified by the low initial value and therefore did not indicate that extension produced the greatest absolute mechanical demand. Flexion consistently generated the highest nucleus pulposus stress, reaching 0.58430 MPa at 794.98 N. In the annulus fibrosus, the motion associated with the maximum stress changed with loading level, and the highest value of 2.66130 MPa occurred during left lateral bending at 794.98 N. These results show that axial load magnitude and motion direction jointly determined the stress response of the two disc components.

Because an intact control model was not included, the present analysis cannot isolate the biomechanical effect of facet resection itself. The findings instead characterize load- and motion-dependent stress responses within the fixed post-foraminoplasty geometry.

4.2. Comparison with Previous Finite Element Studies

The present findings showed both agreement and differences with previous finite element studies. A previous BMI-related finite element study [17] reported that increasing BMI-related loading increased nucleus pulposus stress and deformation after PETD, with flexion producing the most pronounced response. This agrees with the present finding that flexion consistently generated the highest absolute nucleus pulposus stress. Previous foraminoplasty finite element models [16, 24] also demonstrated asymmetric and motion-dependent changes in surgical-level disc stress following superior articular process resection. However, in the present model, the motion producing the maximum annulus fibrosus stress varied across loading levels rather than remaining dominated by one motion direction. These differences may be related to variations in resection geometry, the inclusion or exclusion of annular and nuclear tissue removal, model extent, material definitions, and loading conditions.

Greater compressive loading can alter the mechanical environment of the lumbar segment and increase the demand placed on the intervertebral disc. Previous biomechanical studies have similarly shown that changes in physiological loading influence stress and strain within lumbar structures and affect the mechanical response of the annulus fibrosus [25, 26]. In the present post-foraminoplasty model, the nucleus pulposus was most sensitive during flexion, whereas annular stress showed a stronger dependence on motion direction.

Together, these comparisons indicate agreement in the load-dependent nucleus pulposus response while highlighting model-specific differences in the motion associated with peak annulus fibrosus stress.

4.3. Mechanical Significance of Regional Stress Distribution

The posterior concentration observed in the nucleus pulposus during flexion and the anterolateral and posterolateral concentrations in the annulus fibrosus are mechanically relevant because these regions are commonly associated with disc herniation in experimental and clinical literature [27-29]. In the present model, however, the contours identify regions of elevated stress rather than direct evidence of tissue damage or recurrent herniation.

Based on (Fig. 7), the local stress concentration under flexion loading was predominantly distributed in the posterior region of the nucleus pulposus. This finding may initially appear counterintuitive, as flexion involves forward bending of the spine and the highest von Mises stress was observed during this movement. However, this phenomenon can be explained by the mechanical behavior of the intervertebral disc during flexion, in which compressive stress develops in the anterior annulus while tensile and shear stresses occur in the posterior region. This posterior stress concentration indicates backward migration of the nucleus pulposus during flexion, which increases the loading on the posterior annulus fibrosus. These observations are consistent with previous studies reporting that flexion induces posterior displacement of the nucleus and stress concentration in the posterior annulus, corresponding to the common site of disc herniation [30].

Based on (Fig. 8), the von Mises stress distribution in the annulus fibrosus exhibited two principal concentration areas in the anterolateral and posterolateral regions. This pattern differed from the predominantly posterior concentration observed in the nucleus pulposus during flexion. The broader annular distribution is consistent with load transfer through the circumferential and radial structure of the annulus, while the relative prominence of each region varied with motion direction and axial loading. The posterolateral concentration therefore indicates locally elevated mechanical demand but does not establish annular failure or recurrent herniation.

The maximum annulus fibrosus von Mises stress obtained in the present study was 2.66130 MPa. Experimental tensile testing has reported estimated in situ vertical strengths of approximately 3.9 MPa for the anterior annulus and 8.6 MPa for the posterior annulus [31]. Thus, the maximum stress obtained in the present simulations was numerically below these reported experimental strengths.

Direct equivalence cannot be established because the experimental values represent direction-specific tensile strength, whereas the present study evaluated equivalent von Mises stress using an isotropic linear-elastic annulus model. Annulus fibrosus failure is anisotropic and multiaxial and is strongly influenced by fiber orientation and loading direction [32]. A formal safety factor was therefore not calculated. The reported stress values indicate increased mechanical demand but do not independently demonstrate annular damage or disc failure

4.4. Clinical Implications

The observed load-dependent stress response provides a mechanical context for clinical studies that have associated elevated BMI with less favorable outcomes or recurrent herniation after lumbar surgery [33-35]. The findings support considering body mass together with postoperative movement and loading conditions, while the present simulations do not directly demonstrate annular failure, nucleus extrusion, or recurrent herniation.

The lower equivalent stress observed during extension provides biomechanical context for postoperative movement planning, but it does not establish extension as universally safe or justify a single rehabilitation strategy [36–39]. Rehabilitation should remain individualized and may incorporate core stabilization and proprioceptive training according to symptoms, surgical extent, functional capacity, occupation, and other patient-specific factors [40].

4.5. Study Limitations

This study has several limitations. First, obesity was represented solely through BMI-derived axial compressive loading, while the underlying geometry, material properties, contact definitions, and applied moments were held constant. This controlled approach enabled the load-related mechanical effect of increasing body mass to be isolated; however, it did not capture other obesity-associated biomechanical factors, including altered sagittal alignment, posture, muscle recruitment, muscle-generated forces, fat distribution, and tissue quality. In addition, the simulations were limited to axial compression combined with pure moments and therefore did not reproduce dynamic loading, combined multiplanar movements, or the full complexity of physiological spinal loading. Accordingly, the findings should be interpreted as comparative axial-loading responses within a fixed post-foraminoplasty geometry rather than as comprehensive patient-specific simulations of obesity-related spinal biomechanics.

Second, the study evaluated a single post-foraminoplasty L4–L5 geometry without an intact control model. A single CT-derived anatomy was intentionally used to maintain the same anatomical geometry, material properties, contact definitions, surgical configuration, and boundary conditions across all BMI-derived loading scenarios, thereby allowing the mechanical effect of axial load magnitude to be evaluated under controlled conditions. Similar subject-specific CT-based finite element approaches have been used in previous lumbar biomechanical studies [18, 20]. However, the use of a single anatomy does not capture inter-individual variations in disc dimensions, facet orientation, vertebral morphology, sagittal alignment, and local tissue properties, which may influence both the magnitude and location of peak stress. In addition, because an intact control model was not included, the present analysis cannot independently quantify the biomechanical effect of foraminoplasty relative to the intact lumbar segment. Therefore, the findings are most applicable to the direction and relative pattern of the load-dependent response within the reconstructed anatomy rather than to population-wide absolute stress values. Future studies should repeat the loading protocol across multiple CT-derived geometries and include an intact control model to determine the robustness of the observed patterns.

Third, the surgical model did not incorporate several discectomy-related changes, including annulus fibrosus defects, reduction in nucleus pulposus volume, endplate injury, or other procedure-associated tissue alterations. Several surgical steps were also not explicitly simulated. Consequently, the model represents a simplified post-foraminoplasty condition rather than the complete anatomical changes produced by PTED. No anisotropic tissue-damage or failure criterion was assigned to the annulus fibrosus; therefore, the model could not determine the onset of fiber damage, delamination, or rupture.

Fourth, the finite element model did not incorporate the biomechanical effects of neural decompression. Although nerve decompression is clinically important, its influence on postoperative movement and patient behavior cannot be directly represented using the present structural model. The remaining differences from experimental ROM may also reflect subject-specific geometry, simplified isotropic material behavior, ligament representation, contact definitions, and inter-specimen variability. Future studies should combine multiple subject-specific postoperative anatomies, muscle loading, sagittal alignment, an intact control model, and more comprehensive surgical and damage formulations to improve the physiological relevance of the simulations.