1 Introduction

Degenerative spondylolisthesis (DS) is a condition with varying degrees of slippage between vertebrae caused by degenerative changes in the intervertebral discs and facet joints. DS is a common condition in elderly individuals and one of the most frequently occurring diseases requiring surgical treatment [1]. Degenerative spine pathology in the United States causes a large economic burden, with direct treatment costs estimated at 80 to 100 billion US dollars per year and additional significant indirect costs of disability and lost work productivity [2, 3]. The US National Health Expenditure is projected to increase 1% faster than the gross domestic product per year for 2017 to 2026 [4]. The prevalence of DS in women is approximately six times greater than that in men, and DS is more common in individuals aged > 45 years [5, 6]. Studies have shown that the overall incidence of DS among middle-aged and elderly people in China is 11.20% [7]. The incidence of slippage involving one vertebral body (Meyerding Grade I) is a common degenerative disc disease, accounting for approximately 90% of DS patients [8]. The most common slippage site is L4-L5, accounting for 73% of DS patients, followed by L5-S1 [9, 10]. As our understanding of the pathophysiology and prognosis of DS has improved, it has become increasingly apparent that this condition encompasses a wide spectrum of biomechanical behaviours.

Despite DS being a representative degenerative disc disease, there are relatively few related biomechanical studies. Ramakrishna et al. [11] developed a normal model of the L1-S1 spine and predicted the tendency of vertebral forwards slippage by reducing the disc height. Yaroub et al. [12] manually decreased the disc height and nucleus area and modelled disc height collapse as a 50% reduction in disc height measured in the mid-sagittal plane by moving the L4 inferiorly while maintaining L4–L5 facet articulation. These studies lack objectivity regarding the disease in clinical patients in their model construction, as DS may also involve other degeneration-related structural defects, such as facet joint hypertrophy, intervertebral disc height loss, and adjacent osteophyte formation [13]. At the same time, some studies have used relatively simplistic approaches for the establishment of biomechanical models involving only the adjacent lumbar vertebrae and not incorporating the complex structure of the pelvis beneath the lumbar spine [14]. Changes in the human lumbar spine are not singular or isolated; there is a complex biomechanical relationship between the lumbar spine and pelvis. The onset of DS is also closely related to the lumbar-pelvic position, thus necessitating a more comprehensive organizational structure in simulation studies. Two clinical studies that compared the spinopelvic parameters with global alignment between DS patients and healthy individuals revealed that the pelvic tilt angle (PT) and pelvic incidence angle (PI) were greater in DS patients than in control individuals (i.e., healthy individuals), and the latter was significantly associated with the degree of vertebral slippage [15, 16]. Additionally, muscle tissues play an important role in governing spinal movement and maintaining stability, but current simulation studies that utilize modelling and solving for muscles are still in their infancy [17]. For example, Zhu et al. [18] investigated the effect of back muscle weakness on DS. However, the back muscles were simplified as having a 200 N compressive load, resulting in the incomplete establishment of muscles.

Currently, related clinical studies are unable to quantitatively assess the biomechanical effects of abnormal structures in DS on biological tissues. Understanding the stress‒strain response of biological tissues can help us make more accurate judgements in disease diagnosis and lead to novel therapeutic targets and comprehensive rehabilitation programs, thereby collectively enhancing therapeutic outcomes for patients. The prime example is the lack of available biomechanical evidence to support and guide the design and optimization of specific exercise prescriptions for patients with DS, since backwards bending is a typical risk factor [19]. Therefore, in this study, we established finite element (FE) models of the lumbar-pelvic complex in both DS patients and normal individuals and conducted a comprehensive simulation study from L1-L5 to the entire pelvis to quantitatively assess the biomechanical response differences between the normal lumbar-pelvic composite model and the DS model. The specific objectives of this study are as follows: (1) to establish finite element models representing the lumbar-pelvic complex in both DS patients and healthy individuals, enabling a direct comparison of biomechanical behaviours, and (2) to quantitatively assess the biomechanical response differences, including load distribution, stress‒strain patterns, and joint kinematics, between the normal and DS models to identify critical biomechanical alterations associated with degenerative spondylolisthesis.

2 Materials and methods

2.1 Development of the FE spinopelvic complex model

2.1.1 Geometry reconstruction

The DS spine model and the normal spine model were established based on female patients and asymptomatic adults, respectively, who were recruited for the experiment. The age, height, weight, and demographic data are detailed in Table 1; we recognize that anatomical variations may impact the results; therefore, it is selected patients and asymptomatic adults of the same ethnicity (Asian), gender (female), and with age, height, and weight differences within 10% for finite element modelling to minimize the impact of anatomical variations, which existing researches have reported that these demographic indicators influenced the mechanical responses of a spine motion segment [20,21,22]. This study was approved by the Medical Ethics Committee of the Fourth Clinical Medical College of Guangzhou University of Chinese Medicine (no. K2022-108–02), and informed consent was obtained from the included subjects. The spinopelvic parameters of the DS patients and asymptomatic adults were measured and evaluated by the same physician using X-ray sagittal views, as shown in Table 2. The PI, PT, sacral slope (SS), and lumbar lordosis (LL) angles were all within the normal range for the healthy adults, while the PI, SS, and LL angles of DS patients were all outside the normal range for healthy adults [23,24,25]. The geometric shapes of the bone structures were extracted from the computed tomography (CT) imaging data of the two groups with an axial interval of 0.60 mm and a sagittal and coronal interval of 0.74 mm. The 3D geometry structure was constructed by using image processing software Mimics (version 19.0), which transformed the dicom format image into a digital model. The model was smoothed, amended, and spherized with Geomagic Studio (version 2015). Then, bone structures were used to generate the solid model in SOLIDWORKS CAD software (version 2021). Soft tissue structures such as intervertebral discs, zygapophyseal cartilage, pubic symphysis, and sacrococcygeal joint cartilage were reconstructed based on anatomical studies [26,27,28,29].

Table 1 Demographic features between asymptomatic adult and DS patient
Full size table
Table 2 Comparison of spinopelvic parameters between asymptomatic adult and DS patient
Full size table

2.1.2 Mesh generation and modelling

Mesh generation and modelling were performed with HyperMesh software. The generation process in combination with an appropriate mesh resolution allows preservation of morphological features of the biological tissue in the model. The overall FE spinopelvic complex model was developed using Ls-dyna codes (version 11.1.0 × 64, simple precision, SMP), as shown in Fig. 1. The element type selections are listed in Table 3, and numbers of element (nodes) between two models are listed in Table 4.

Fig. 1
figure 1

a View of the normal spine in the 3D reconstructed model. b Spinal function unit of the lumbar spine. c View of the DS spine in the 3D reconstructed model. d L4/5 segment with degenerative change (slippage forward) visible

Full size image
Table 3 Summary of material and element properties used for the spinopelvic complex FE model
Full size table
Table 4 Numbers of element (nodes) between normal spinopelvic complex and DS spinopelvic complex models
Full size table

Functional spinal unit (FSU). Bony structures were modelled by combining shell and solid elements. The cortical bones of each vertebra were modelled with triangular elements to accurately fit the anatomical features, while their thicknesses were based on cadaver data summarized in the literature [30]. The trabecular bone was modelled as tetrahedral elements. The zygapophyseal joint cartilages and bony endplates of the model were generated based on vertebral landmarks and anatomical descriptions [31, 32]. Cappetti et al. [33] reported that geometric parameters significantly affect disc mechanical behaviour, such as sagittal thickness (\({s}_{d}\)), width (\({\overline{w} }_{d}\)), and height (\(H\)) of the disc. In our study, the intervertebral disc modelling parameters are sagittal thickness (\({s}_{d}\)) of 28.79 ± 1.03 mm, width \({\overline{w} }_{d}\) of 43.63 ± 2.24 mm, and height (\(H\)) of 7.96 ± 1.80 mm, which are comparable in geometric size to the spinal segment model parameters reported by Cappetti et al. [33] with differences all within 5%. The facet joints were modelled with triangular elements to accurately fit the anatomy features. The frictionless contacts were defined between the superior and inferior articular cartilages, which were tied to the lateral masses of the vertebrae. The intervertebral discs were modelled as composites of hexahedral bulk element and quadrilateral membrane elements, tied to the vertebral endplates. The intervertebral discs included the components of annulus fibrosus matrix, nucleus pulposus, and annular fibres, which were respectively modelled with hexahedral bulk elements or quadrilateral membrane elements [34]. The intervertebral discs were tied to the adjacent vertebral endplates. The central nucleus pulposus, the anulus fibrosus was four elements in the cranial–caudal direction and two elements in the radial direction. On each layer of the anulus in the radial direction, two layers of membrane elements were superimposed, to represent the nonlinear and orthotropic characteristics of the anulus fibrosus [35, 36]. The cartilaginous endplates were attached to the disc surfaces in the cranial-caudal direction.

Ligaments

In the overall lumbar spine (LS), from L1 to L5, the anterior longitudinal ligament (ALL), posterior longitudinal ligament (PLL), ligamentum flavum (LF), interspinous ligament (ISL), and supraspinal ligament (SSL) were modelled with quadrilateral membrane elements based on anatomical descriptions [37,38,39]. Facet capsular ligaments (FCLs) were modelled with 1D spring elements, and facet cartilage was defined by shell elements. The modelling details of the hip joint capsule were described by Anantha-Krishnan et al. [40]; the joint capsule was modelled as a cylindrical sleeve that originated from an ellipse around the acetabular rim and inserted into an ellipse around the femoral intertrochanteric line. Then, we generated a series of nodes based on anatomical sites [41]. The corresponding nodes in adjacent lines were connected to form 4-node quadrilateral elements along the length of the capsule. The hip joint capsules were subdivided into four longitudinal sectors to approximate the anatomical details.

Muscles

As shown in Fig. 2, the spinopelvic muscles were developed and connected to adjacent bony structures based on a previously published study [36]. Seventy-four of Hill-type 1D elements were used to represent 6 muscle pairs, including the longissimus thoracis pars thoracis, iliacus muscle, quadratus lumborum, musculus iliocostalis, psoas major, and psoas minor, in the spinopelvic complex. The muscle fascicles were redirected at the intermediate vertebrae to better follow the curvature of the spinopelvic structure and upper femurs.

Fig. 2
figure 2

View of the spinopelvic muscles in the FE model (a ventral muscles and b lateral and dorsal muscles)

Full size image

2.1.3 Material properties

The material properties of the model were based on previous studies, as shown in Table 3. Bone material parameters referenced from classical FE modelling studies, including the whole lumbar spine models by Park et al. [42], Little et al. [43], Rohlmann et al. [44], Goel et al. [45, 46], the ViVA human body model [47, 48], and the GHBMC (Global Human Body Models Consortium) model [49, 50]. Isotropic and kinematic hardening plasticity was used to model the cortical bones, cancellous bones, and bony endplates of the lumbar vertebrae, pelvis, and femurs [36, 51]. Referring to the studies by Zhang et al. [52], the elastic modulus of cortical and trabecular bone was reduced by 33% to represent the degenerative changes in bone mineral density. The nucleus pulposus was modelled with hypoelasticity material [53]. The annulus was built of a multipart material consisting of an annulus fibrosus matrix in a hill foam material [54, 55] and annular fibres with nonlinear orthotropic characteristics [56]. The spinopelvic ligaments were represented by different anatomical thicknesses using orthotropic nonlinear elastic material, as indicated by previously published studies [49,50,51]. Hip joint capsules and symphysis pubis were modelled by linear elastic materials with elastic moduli of 50 MPa and 70 MPa [57,58,59], respectively. The facet cartilage was also modelled using a linear elastic material with an elastic modulus of 10 MPa [60]. Spinopelvic muscles were added to the spinopelvic complex FE model and defined by the improved Hill model. Detailed data such as the origin and insertion points, physiological cross-sectional areas, and muscle segment parameters are available in previous comprehensive dissection studies [36, 61]. Muscle passive properties were determined according to previous studies [48, 62].

2.2 Model validation

The validation matrix of the spinopelvic complex model is illustrated in Table 5. The model was first verified via anatomical measurements. Then, it was validated step by step from the isolated tissues to the segments. For all the conducted simulations, gravity loading (a volume load of 9.81 m/s2) was included. The simulation details are described as follows.

Table 5 Summary of validation matrix for the spinopelvic complex FE model
Full size table

2.2.1 Compression test

First, as described by Wagnac et al. [63, 64], a uniform axial displacement rate of 0.1 m/s (low dynamic loading) and 1.0 m/s (high dynamic loading) for the superior nodes of the vertebral body was used to evaluate the stiffness of the vertebral body. The simulated force–displacement curve was compared with that of the experimental corridor.

2.2.2 Tension tests

Tension tests were performed on the ALL, PLL, LF, ISL, SSL, and FCL. The respective force‒displacement curves were compared with the experimental data from Wagnac et al. [63, 64].

2.2.3 6-DOF moment of functional spinal units

Next, the motion characteristics of the model were verified in three rotation directions as described in previous studies [44, 65]. For the entire lumbar spine and each FSU, moments were applied to the superior portion of the L1 at 3.75 Nm and 7.50 Nm for flexion–extension (F-E), lateral bending (LB), and axial rotation (AR), as described in previous studies [44, 65]. The respective moment–rotation simulation curves were compared with the experimental ranges.

2.2.4 Shear test for the DS model

Lateral shear tests were performed on the cranial vertebra relative to the caudal vertebra of L4, since L4/5 is the slippage segment in the DS model. The pure LB moments were applied at 5 Nm in accordance with the load increments in the study by Melnyk et al. [66]. In addition, the L4/5 DS model was subjected to anterior shear and posterior shear to the cranial vertebra. The moment–rotation simulation curves and force–translation curves were compared with the experimental data [66].

All simulations were run with single-precision LS-DYNA codes (version 11.1.0 × 64, simple precision, SMP). The bilateral femurs were constrained in all directions. General contact was enforced for the entire model. The motion characteristics of the segments and the intervertebral discs as well as the vertebral loads were analysed. A global coordinate system with X to the left, Y backwards, and Z upwards was utilized. Finally, to explore the biomechanical differences in DS involvement under different loading modes, the responses of the spinopelvic kinematics and important adjacent tissues or segments were compared between the normal spinopelvic FE model and DS spinopelvic model.

3 Results

3.1 Validation of the spinopelvic FE model

3.1.1 Compression test

The vertebrae were subjected to compressive tests, simulating force‒displacement curves that ranged from high-dynamic (HD) to low-dynamic (LD) conditions. The results of this comprehensive evaluation are illustrated in Fig. 3a, b. Under HD conditions, a substantial reaction force of 5605.3 N was recorded at a displacement of 3.7 mm. The simulated force‒displacement curves and the experimental values were similar, exhibiting a monotonic increasing trend, with the peak values being nearly consistent. Under LD conditions, a lower reaction force of 1902.0 N was observed at a displacement of 6.5 mm. The simulated force‒displacement curve for LD was scattered within the experimental range.

Fig. 3
figure 3

a High dynamic compression validation of the normal spinopelvic model compared with Wagnac et al. [63, 64]. b Low dynamic compression validation of the normal spinopelvic model compared with Wagnac et al. [63, 64]

Full size image

3.1.2 Isolated tissues

In spinopelvic ligament tension tests, most of the force–displacement curves for the normal model were within the experimental corridor concerning both curve shapes and peak values, as shown in Fig. 4. Most simulated force‒displacement curves could be described as a distinctive J-shaped curvature up to the peak load and exhibited regions of toe, linear, and traumatic motion in the process of tension. The PLL simulation curve exhibited an asymmetric sigmoidal shape up to the peak load. The validated results delineated the range of the ultimate load capacity of the spinopelvic ligaments, which spanned from approximately 59.8 to 411.5 N. This wide range reflects the significant variability in the mechanical strength of different ligaments within the spinopelvic framework, highlighting the complex interplay of forces that these structures can withstand in vivo.

Fig. 4
figure 4

a ALL tension validation of the normal spinopelvic model compared with Wagnac et al. [63, 64]. b PLL tension validation. c LF tension validation. d FCL tension validation. e SSL tension validation. f ISL tension validation

Full size image

3.1.3 6-DOF moment of functional spinal units

For entire spinopelvic FSU rotation analysis, this study employed a comprehensive 6-DOF rotational analysis to elucidate the biomechanical responses to pure moments. The corresponding full moment–rotation responses, encapsulated in the normal spinopelvic model, are meticulously documented in Figs. 5 and 6, including the F-E, LB, and AR loading conditions. The pure moments were set as 3.75 Nm and 7.5 Nm for assessing the ROM within the spinal segments, providing a detailed comparative analysis against established experimental corridors. In the F-E, LB, and AR loading conditions with pure moments of 3.75 Nm and 7.5 Nm, most of the ROMs were within the experimental corridors for both curve shapes and peak values. This congruence underscores the reliability and accuracy of the simulation model in replicating the physiological and biomechanical properties of spinopelvic FSUs under dynamic loading conditions. Additionally, the ROMs of the separated FSUs in the entire spinopelvic FE model were all within one standard deviation of the experimental data [65] for loading with pure moments of 7.5 Nm, as shown in Fig. 7.

Fig. 5
figure 5

Moment–displacement curve validation of the normal spinopelvic model for loading with pure moments of 3.75 Nm in a flexion–extension, b lateral bending, and c axial rotation compared with Rohlmann et al. [44]

Full size image
Fig. 6
figure 6

Moment–displacement curve validation of the normal spinopelvic model for loading with pure moments of 7.5 Nm in a flexion–extension, b lateral bending, and c axial rotation compared with Rohlmann et al. [44]

Full size image
Fig. 7
figure 7

7.5 Nm moment histogram validation of each segment (L1-S1) in the normal spinopelvic model under a flexion–extension, b lateral bending, and c axial rotation compared with Cook et al. [65]

Full size image

3.1.4 Shear test for the DS model

In the intricate field of biomechanical engineering, particularly in the analysis of the DS model, the shear test has emerged as a critical methodological approach to assess mechanical behaviour under varying load conditions. This study compared moment–rotation simulation curves and force–displacement curves against experimental data [66], thereby providing a robust model for evaluating the predictive accuracy and reliability of the DS spine under shear stress conditions. The comparative analyses are graphically presented in Fig. 8. The moment–rotation simulation curve (Fig. 8a) exhibited an asymmetric S-shape, and the majority of the curve fell within the experimental corridor. The anterior–posterior translation had a relatively linear relationship with the shear force across the shear loading conditions, as shown in Fig. 8b. Although the peak value of the anterior shear curve in the DS model was slightly lower than the experimental peak value, the tendency was consistent with that of the experimental curve. The simulation results and experimental findings aligned closely, particularly in the context of posterior shear loading.

Fig. 8
figure 8

a Lateral shear validation of the DS model compared with an in vitro model of DS from Melnyk et al. [66]. b Anterior–posterior shear validation of the DS model compared with the experiment data from Melnyk et al. [66]

Full size image

3.2 Comparison of the normal spine model and DS spine model

The peak stress of the lumbar isthmic-cortical bones, intervertebral discs, and facet joints; the peak strain of ligaments; the peak force of muscles; and the percentage difference in the ROM during F-E, LB, and AR motions were analysed and compared between the normal spine model and the DS spine model.

3.2.1 Lumbosacral intervertebral disc stress

The peak stress of the intervertebral discs under F-E, LB, and AR conditions was greater than that in the normal model. As shown in Fig. 9a, in the comparison of the peak stress on the L4/5 intervertebral disc (Fig. 9a), under the influence of L4/5 segment slippage, the stress in the DS spine model was greater than that in the normal model; under F-E, LB, and AR conditions, the stress was 1.6, 1.1, and 6.6 times greater than that in the normal model, respectively. In the normal spine model, the least stress on the intervertebral disc occurred under the AR condition, whereas in the DS spine model, the highest stress was observed under the AR condition. Following disc slippage, there is a decrease in the segment’s ability to resist AR motion, a pattern that also becomes apparent during F-E motion. A comparison of the peak stress in the L5/S1 intervertebral disc between the DS and normal spine models revealed minor differences, as shown in Fig. 9b. However, under LB and AR conditions, the peak stress of the normal model at the L5/S1 disc was still 16.18% and 14.84% lower, respectively, than that of the DS spine model.

Fig. 9
figure 9

Comparison of peak stress in a L4/5 intervertebral disc and b L5/S1 intervertebral disc between normal spine model and DS spine model

Full size image

3.2.2 L4 cortical bone von Mises stress and typical stress distribution nephograms

As illustrated in Fig. 10, we compared the peak stress of the L4 cortical bone between the two models under F-E, LB, and AR conditions. In terms of the distribution of peak stress in the L4 cortical bone under different conditions, the values for the DS spine model were consistently greater than those for the normal model; the peak stress was 1.2, 2.1, and 1.7 times greater under F-E, LB, and AR conditions, respectively. The peak stress presented by the L4 cortical bone in the DS spine model was the greatest during axial rotation. According to the typical stress distribution nephograms of flexion motion (Fig. 11), the overall stress region of the DS spine model was greater than that of the normal spine model. The former showed more high-stress areas than did the latter, especially in anatomical structures prone to stress concentration, such as the facet joint isthmus, where the DS spine model exhibited greater peak stress than did the normal spine model. Additionally, forwards slippage of the L4/5 pathological segment caused a shift in the high-stress distribution nephogram of the vertebrae towards the vertebral pedicle and anterior edge of the lamina.

Fig. 10
figure 10

Comparison of peak stress in L4 cortical bone between DS spine model and normal spine model

Full size image
Fig. 11
figure 11

Comparison of typical von Mises stress (GPa) distribution nephograms (flexion motion) between two models. a Top view of L4 in normal spine model. b Top view of L4 in DS spine model. c Rear view of L4 in normal spine model. d Rear view of L4 in DS spine model

Full size image

3.2.3 Percentage differences in L4/5 and L5/S1 ROMs

Comparison of the percentage differences in mobility levels of the L4/5 and L5/S1 segments between the normal spine model and the DS spine model clearly illustrates the disparities in segment stability between the models. As shown in Fig. 12, both the L4/5 and L5/S1 segments were significantly greater in the DS spine model than in the normal spine model, indicating that slippage pathology affects segment stability. The greatest impacts primarily occurred during AR and LB motions, followed by the F-E motion. This finding suggested that after anterior vertebral slippage, the pathological segment mainly exhibited reduced abilities to resist AR and LB motions. The comparison between the two models in terms of the percentage difference in the L5/S1 ROM reflected a similar trend, although segments without DS did not show compromised stability during F-E motion.

Fig. 12
figure 12

Comparison of percentage difference of L4/5 and L5/S1 ROM between the normal spine model and DS spine model

Full size image

3.2.4 Peak stress of the L4/5 facet joint and peak strain of the spinal ligaments

Given the notable differences in biomechanical responses to LB and AR motions between the normal spine model and DS spine model and considering the vital anatomical function of the facet joints in resisting LB and AR motions, we compared the stress differences in the L4/5 facet joints of the two models. As shown in Fig. 13, a substantial disparity in facet joint stress was observed between the two models, with the DS spine model consistently showing greater stress than the normal model across the three conditions, specifically, 1.3, 6.6, and 1.9 times greater for the F-E, LB, and AR loading conditions, respectively.

Fig. 13
figure 13

Comparison of L4/5 facet joint peak stress between normal spine model and DS spine model

Full size image

Compared to the loading conditions of F-E, LB, and AR (Fig. 14), the peak strain values of the ligaments in the DS spine model were greater than those in the normal spine model. The peak strain of the PLL in the DS spine model was the highest among all the tested spinal ligaments and was notably greater than that in the normal spine model. The PLL peak strain of the normal model was only 1/10 that of the DS model, indicating that the PLL is the most susceptible ligament that is significantly affected after anterior vertebral slippage. In addition, a pronounced disparity was also evident in the peak strain of the ISL between the two models. A notable difference in the peak strain of the ISL was observed between the two models, particularly for the F-E and AR motions.

Fig. 14
figure 14

Comparison of L4/5 ligaments peak strain between the normal spine model and DS spine model. (a) Under FE loading condition. (b) Under LB loading condition. (c) Under AR loading condition

Full size image

3.2.5 Muscle tensile force in the lumbosacral region

As shown in Fig. 15, a comparison of forces under six conditions revealed that the trends in muscle tensile forces were generally similar between the two models. However, most muscle tensile forces in the lumbosacral region of the DS spine model were greater than those in the normal model, except the iliacus muscle. The psoas muscle and iliocostalis muscle exhibited the greatest tensile forces in the two spinopelvic complex models. During spine extension, the psoas muscle in the DS spine model was subjected to 23.2% greater tensile force than that in the normal spine model. Under LB loading conditions, the iliocostalis muscle in the DS spine model also exhibited a 42.2% greater tensile force than that of the normal spine model. Compared with those of the normal spine model, the tensile force of the quadratus muscle of the DS model increased by 34.7%. Overall, the DS model exhibited greater muscle tensile forces in the lumbosacral muscle groups during F-E and LB motions than did the normal model, indicating that L4 anterior slippage and changes in lumbosacral-pelvic alignment significantly affected the biomechanical response of the muscles. The comparative muscle nephograms in Fig. 16 illustrate the tensile forces on the psoas during extension in both models, revealing that the muscle tension in the DS spine model exceeded that in the normal model, with a force concentration on the DS pathological segment and its proximal segments.

Fig. 15
figure 15

Comparison of muscle (iliacus, quadratus, longissimus thoracis, psoas, and iliocostalis) tensile force of lumbosacral region between normal spine model and DS spine model

Full size image
Fig. 16
figure 16

Comparative muscle nephograms of the muscle tensile forces on the psoas during extension between a normal spine model and b DS spine model

Full size image

4 Discussion

Based on the development and validation of a new spinopelvic complex FE model with muscle modelling, the present study simulated an in vivo response of the human spinopelvic system and analysed the effects of DS with abnormal spinopelvic parameters on lumbosacral kinematics. Due to the pathological changes involved in vertebral slippage in the DS model and the corresponding changes in spinopelvic positional parameters, alterations were observed in muscular and skeletal structures, tissue strain, stress levels, segmental ROM, and load concentration locations. Compared to the results of cadaveric experiments, the validated analysis demonstrated the biofidelity of our normal spine model, as shown in Figs. 3, 4, 5, 6, and 7. Moreover, limited studies have recorded scattered data, which are also well correlated with the present DS model. As shown in Fig. 8, the DS model showed segmental moment–rotation curves and force–displacement curves mostly within the experimental corridor and consistent with the curve trend.

While several studies have addressed the biomechanical effects of DS, previous computational models have been limited to isolated lumbar spine studies and simplified musculature. Yaroub et al. [12] reconstructed the FSU of the fourth and fifth lumbar vertebrae with spondylolysis and reported that AR motion was associated with the most significant increase in the ROM in an isthmic spondylolisthesis model, which was consistent with our finding that the DS spine model displayed hypermobility in terms of flexion and rotation. Liu et al. [67] also developed a L3-S1 FE model and demonstrated that intradiscal pressure was increased in the DS spine model, indicating a greater possibility of disc degeneration. However, they did not use a suitable method of validation for the spondylolisthesis model, as it was developed by extending the isthmus of the normal lumbar spine, which may not accurately reflect the morphological characteristics of lower-grader DS. Only one study evaluated the effects of back muscles (simplified as a 200 N compressive load) on DS, which could impact the robustness of the conclusions [18]. A more detailed model including additional muscles would be beneficial for determining the biomechanical effects of abnormal DS structures on biological tissues. In this study, a DS spine model was reconstructed from a real DS patient whose age, height, and weight were within 10% of those in the normal spine model. However, significant differences in the spinal-pelvic parameters existed between the two models, with the PI, SS, LL, and PT parameters in the DS spine being greater than those in the normal spine, while the first three parameters even exceeded the normal ranges. Changes in the spinopelvic structure of DS cause alterations in adjacent skeletal muscles, combined with progressive degenerative changes in the pathological segments. Additionally, changes in the line of force applied to the spinopelvic structure, which causes changes in the sagittal plane alignment of the spine, have been implicated in the development and progression of DS. Barrey et al. [68] proposed that an increase in the SS and PI angles from the normal reference values was predictive of DS. In agreement with these findings, Curylo et al. [69] reported a significantly greater PI among fifty-three patients with DS than among a control group. Labelle et al. [70] also reported greater LL values among patients with DS, with an increase in this parameter being positively correlated with the grade of spondylolisthesis. These studies confirm a relationship between DS and spinopelvic parameters of sagittal plane alignment, with an increase in the PI and LL being specifically associated with DS. Therefore, when numerically modelling and analysing the DS model, it is crucial to consider the positional parameters of the lumbar pelvis. Global sagittal balance has unique biomechanical effects on the biological tissues of pathological segments and their adjacent segments in the DS model.

Regarding the intervertebral discs, the two models showed substantial differences concerning both peak stress values, as shown in Fig. 9. Under the influence of L4 forwards slippage, the peak stress of the L4/5 intervertebral disc in the DS spine model under the AR condition was 6.6 times greater than that in the normal spine model. Following disc slippage, the ability of the segment to resist AR motion decreases, and this tendency also becomes apparent with slippage at the L5/S1 segment. Following degenerative changes in the intervertebral disc, the height of the disc space decreases, the structure of the endplate deforms, and the ability of the disc to absorb shocks and distribute loads evenly across the endplate decreases, ultimately reducing the resistance of the disc to AR motion. The primary biomechanical function of the intervertebral disc is to resist compressive forces, and its overall structure is not conducive to withstanding other types of loads; it is especially poor at tolerating torsional stress [71]. Additionally, the increased peak stress in the intervertebral discs of spondylolisthesis patients and the adjacent segments also influences posterior facet joints, which are crucial anatomical structures in the functional spinal unit opposing LB and AR motions. As DS often accompanies an increase in lumbar pelvic parameters, the increased shear forces in the segment also increase the stress on facet joints and intervertebral discs. Labelle et al. [72] reported that a linear association existed between the PI and spinal spondylolisthesis. Lv et al. [73] reported that patients with increased PI values are at a greater risk of DS, and an increase in PI leads to elevated mechanical stress on lumbar facet joints. Our results also indicated that the facet joint stress in the L4/L5 pathological segments of the DS spine model was notably greater than that in the normal spine model.

Increased stress on the facet joints and intervertebral discs suggested a greater ROM during LB and AR motions in the DS spine model, as shown in Fig. 12. The reported data indicated that segmental motion was greater in patients with moderate facet joint and disc degeneration than in normal subjects [74]. Increased mechanical stress on the three-joint complex in the functional spinal unit leading to segmental instability is an important cause of low back pain, potentially accelerating tissue degeneration, including hypertrophy of the facet joints and further reduction in disc space height, resulting in the redistribution of stress as the biological tissues attempt to achieve segmental adaptive stabilization [75]. Our findings also suggest that the ROM in the F-E motion of segments following vertebral forwards slippage is only slightly greater than that in the normal spine model. Changes following DS also cause stress alterations, as illustrated by the typical stress distribution nephograms in Fig. 11. Under flexion motion, the DS spine model showed more high-stress areas than did the normal spine model, especially in anatomical structures prone to stress concentration, such as the facet joint isthmus, where the DS spine model exhibited greater peak stresses than did the normal spine model. Additionally, forwards slippage of the vertebra resulted in increased strain on the posterior spinal ligaments, notably the PLL and ISL. Overall, the trend in changes in peak ligament strain indicated that the DS spine model exhibited obviously greater strain than the normal spine model.

The natural history of DS is generally favourable. Only 10–15% of patients seeking treatment will eventually undergo surgery [76]. Therefore, for patients with DS, conservative treatment remains the primary therapeutic approach in current clinical practice. Overuse of the stabilization muscles may cause or aggravate forwards slipping in DS patients. However, little biomechanical evidence is available to explain the mechanism underlying exercise therapy [18]. Merely strengthening the back muscles appears to be insufficient because DS is often accompanied by an anterior pelvic tilt and an increased lordotic angle of the lumbar spine, resulting in the shortening and tightening of anterior muscles such as the psoas. As these weakened muscles stay in that position, they need to work extra hard to stabilize the pelvis and core, causing them to tighten. This study demonstrated that changes in the pelvic tilt affect the length of the muscles around the spine and hip joint. These muscles often become inactive due to the prolonged sitting that occurs in modern life and work. In this study, we found that most muscle tensile forces in the lumbosacral region of the DS spine model were greater than those in the normal model. During spine extension, the psoas in the DS spine model was subjected to 23.2% greater tensile force than that in the normal spine model. Compared with that of the normal spine model, the tensile force of the quadratus muscle of the DS model increased by 34.7%. This result indicates that during the muscle strengthening process in DS patients, it is necessary to strengthen both the ventral and dorsal musculature of the spine-pelvis, whereas most current clinical practices focus primarily on the dorsal muscles [77], inadvertently exacerbating the imbalance in the antagonist muscle groups. In fact, weakness in the abdominal muscles is observed in most patients with low back pain. The anatomical position of the psoas major is very close to the axis of motion of the lumbar spine and hip joint, thereby aiding in maintaining the stability and coaxial alignment of the lumbar spine and hip joints [26, 78]. A weakened psoas in the hip joint can cause problems with balance and posture control during activities. Kim et al. [79] reported that a change in the pelvic tilt influences the alignment of the spine and causes great stress on the lumbar vertebrae, which is consistent with our simulated results. Therefore, this study provides quantitative data to emphasize the importance of exercising the ventral muscle groups in DS patients, aimed at correcting the pelvic position to improve posture.

The limitations of this study should be noted. The absence of active muscular contraction prevented direct evaluation of lumbar spine motion [80, 81]. However, defining active muscle forces in DS patients is particularly challenging due to the complexity and variability of muscle activation patterns in this condition. Neglecting this effect is a simplification used in the present work, which should be considered in future improvements of our model. Additionally, changes in size and shape of bones and tissue may affect magnitude and even quality of the results reported. In this study, by controlling these critical influencing variables, we aim to minimize differences in vertebral size, shape, and other morphological characteristics. The modelling of the DS spine reflects clinical authenticity derived from actual DS patient, and the model has been validated through cadaveric experiments. This methodology obtains model’s biofidelity to a certain degree. In the future, we will utilize spatial interpolation methods, such as Kriging, to establish the FE models. Finally, the vertebrae [82, 83], intervertebral disc, and other tissue materials in the FE model of the lumbar-pelvic complex do not represent the state-of-the-art constitutive equations, which may affect the quantitative data of the simulation results. A notable example includes the annulus fibrosus, which are recently based on biphasic or biphasic-swelling theories [84,85,86,87,88]. These models offer the advantage of simulating fluid–solid interactions and osmotic pressure. The annulus fibrosus modelling of this study is the lack of comprehensive integration among fluid–solid interaction (permeability), the impact of biochemical components (hydration), regional variations in structural arrangement, and mechanical properties of both collagen and elastic fibres (fibre reinforcement). In this study, degenerative spondylolisthesis mainly involves structural disc and facet joints degeneration, while material degeneration, which is also critically important, was not addressed. However, this study is able to show the general trends of load transfers between the DS spine and the healthy spine. More reliable material properties and more accurate geometries by utilizing in vivo experiments, such as Amirouche et al. [89], focusing on DS subjects are needed in order to increase the biomechanical predictive capability and clinical utility of any future the FE models.

5 Conclusion

The present study established, validated, and compared normal spine and DS models based on detailed FE spinopelvic complex models with muscle modelling. The obvious simulation differences indicated greater peak stress in the intervertebral discs of the pathological segment and adjacent segments in DS patients with than in those of the normal spine, as well as overall greater ligament strain and facet joint stress. This is particularly evident during axial rotation of the spondylolisthesis spine. Biomechanical differences originate not only from alterations in the spondylolisthesis segment but also from an overall imbalance in the spinopelvic complex. Forwards vertebra slippage and greater spinopelvic parameters also result in more high-stress areas, such as the isthmus region, and stress-concentrated locations are evidently transferred forwards. Finally, this study provided a quantitative biomechanical comparison to emphasize the importance of exercises targeting the psoas and other ventral muscles for DS patients to improve the pathological state of anterior pelvic tilt and prevent the progression of DS.