Front. Bioeng. Biotechnol.Frontiers in Bioengineering and BiotechnologyFront. Bioeng. Biotechnol.2296-4185Frontiers Media S.A.111075210.3389/fbioe.2023.1110752Bioengineering and BiotechnologyOriginal ResearchAnatomical parameters alter the biomechanical responses of adjacent segments following lumbar fusion surgery: Personalized poroelastic finite element modelling investigationsNikkhoo et al.10.3389/fbioe.2023.1110752NikkhooMohammad123†ChenWen-Chien24†LuMeng-Ling25FuChen-Ju26NiuChi-Chien27LienHen-Yu1ChengChih-Hsiu12*1School of Physical Therapy and Graduate Institute of Rehabilitation Science, College of Medicine, Chang Gung University, Taoyuan, Taiwan2Bone and Joint Research Center, Chang Gung Memorial Hospital, Linkou, Taiwan3Department of Biomedical Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran4Department of Orthopedic Surgery, Chang Gung Memorial Hospital, Taoyuan, Taiwan5Department of Orthopedic Surgery, Chang Gung Memorial Hospital, Kaohsiung, Taiwan6Division of Emergency and Critical Care Radiology, Chang Gung Memorial Hospital, Linkou, Taiwan7Department of Orthopedic Surgery, Chang Gung Memorial Hospital, Linkou, Taiwan
Edited by:Lizhen Wang, Beihang University, China
Reviewed by:Zhongjun Mo, National Research Center for Rehabilitation Technical Aids, China
ChengFei Du, Tianjin University of Technology, China
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Introduction: While the short-term post-operative outcome of lumbar fusion is satisfying for most patients, adjacent segment disease (ASD) can be prevalent in long-term clinical observations. It might be valuable to investigate if inherent geometrical differences among patients can significantly alter the biomechanics of adjacent levels post-surgery. This study aimed to utilize a validated geometrically personalized poroelastic finite element (FE) modeling technique to evaluate the alteration of biomechanical response in adjacent segments post-fusion.
Methods: Thirty patients were categorized for evaluation in this study into two distinct groups [i.e., 1) non-ASD and 2) ASD patients] based on other long-term clinical follow-up investigations. To evaluate the time-dependent responses of the models subjected to cyclic loading, a daily cyclic loading scenario was applied to the FE models. Different rotational movements in different planes were superimposed using a 10 Nm moment after daily loading to compare the rotational motions with those at the beginning of cyclic loading. The biomechanical responses of the lumbosacral FE spine models in both groups were analyzed and compared before and after daily loading.
Results: The achieved comparative errors between the FE results and clinical images were on average below 20% and 25% for pre-op and post-op models, respectively, which confirms the applicability of this predictive algorithm for rough pre-planning estimations. The results showed that the disc height loss and fluid loss were increased for the adjacent discs in post-op models after 16 h of cyclic loading. In addition, significant differences in disc height loss and fluid loss were observed between the patients who were in the non-ASD and ASD groups. Similarly, the increased stress and fiber strain in the annulus fibrosus (AF) was higher in the adjacent level of post-op models. However, the calculated stress and fiber strain values were significantly higher for patients with ASD.
Discussion: Evaluating the biomechanical response of pre-op and post-op modeling in the non-ASD and ASD groups showed that the inherent geometric differences among patients cause significant variations in the estimated mechanical response. In conclusion, the results of the current study highlighted the effect of geometrical parameters (which may refer to the anatomical conditions or the induced modifications regarding surgical techniques) on time-dependent responses of lumbar spine biomechanics.
personalized modelingfinite element analysisposterior lumbar fusionadjacent segment diseasespine biomechanicsNational Science and Technology Council10.13039/501100020950Chang Gung Memorial Hospital, Linkou10.13039/501100005795Chang Gung University10.13039/5011000028361 Introduction
The posterior instrumentation with rigid-rod fusion is considered the gold standard surgical treatment of pathologies such as lumbar instability, spinal stenosis, spondylolysis, and spondylolytic spondylolisthesis (Serhan et al., 2011; Kim and Choi, 2018). Although the short-term post-operative (post-op) outcome of lumbar fusion is satisfying for most patients, adjacent segment disease (ASD) can be prevalent in long-term clinical observations (Kim et al., 2016). This long-term phenomenon mainly affects the adjacent intervertebral discs (IVDs), however, it can reveal instability, retro-spondylolisthesis, and fracture in adjacent vertebrae, as well (Liang et al., 2014; Kim et al., 2016).
Different risk factors have been proposed for increasing the chance of ASD development in patients [such as age, sex, body mass index (BMI), and osteoporosis], non-etheless, the post-op ASD alteration may possibly be a result of induced modifications in lumbosacral spine lordosis angle, kinematics, and kinetics (Helgeson et al., 2013; Vergroesen et al., 2015; Ebrahimkhani et al., 2022). Hence, characterizing the effect of lumbar fusion on biomechanical responses of adjacent levels could be beneficial for surgeons towards improved surgical planning and enhanced clinical outcomes. While various in-vitro experimental (Phillips et al., 2006; Erbulut et al., 2013; Beckmann et al., 2020) and clinical investigations (Kim et al., 2012; Hsieh et al., 2016; Zhang et al., 2018a) have been performed to compare the outcome of using different posterior fusion devices for treatment of the lumbar spine diseases, there is no access to a non-invasive technique to evaluate the post-op alterations in the spinal biomechanics. Finite element (FE) analysis could be similarly employed as a conventional predictive approach for clinical investigations due to its ability to represent the complex systems and predict their response (Dreischarf et al., 2014).
Despite the proven success of FE analyses for investigation of the lumbar spine post surgeries (Zhang et al., 2018b; Zhang et al., 2018c; Zhao et al., 2018), their application in clinical investigations, as an assistive tool, may be questioned by clinicians. The differences in the spine anatomical geometry may possibly cause uncertainty in FE model outputs and limit the reliability of achieved predictions (Laville et al., 2009). Geometrically personalized FE modeling, regardless of its simplifications, can provide a modular tool for clinical studies that include patient-specific characterization using clinical images to account for the variability between different patients (Nikkhoo et al., 2020; Ebrahimkhani et al., 2022).
Hence, to fill the gap of knowledge, it might be valuable to investigate if inherent geometric differences among patients can significantly alter the biomechanics of adjacent levels post-surgery. We recently developed a geometrical personalized FE modeling technique in which the detailed time-dependent fluid-solid interactive response was considered to enhance prediction under both static and dynamic loading conditions (Nikkhoo et al., 2021). The main objective of this study was to utilize this modeling technique to evaluate if the geometrically personalized FE modeling can predict the alterations in lumbar adjacent levels post-fusion surgery. For this purpose, a prospective, non-randomized cohort study was performed in which the patients underwent one-level lumbar interbody fusion and it is hypothesized that the FE modeling technique may identify remarkable variations in adjacent segment kinematics and kinetics between patients without ASD and patients with ASD.
2 Materials and methods2.1 Geometrical personalized FE modeling of the pre-operative lumbosacral spine
The pre-operative geometries of the lumbosacral spine (L1-S1) of 30 patients were generated from lateral and anterior-posterior (AP) X-ray radiographs (Age: 64.8 ± 8.1 years, BMI: 26.2 ± 3.5 kg/m2, 25 females and 5 males) using a previously developed validated geometrical modeling procedure (Nikkhoo et al., 2020) (Figure 1). The X-ray radiographs were selected from a prospective, non-randomized cohort study in which the patients underwent one-level lumbar interbody fusion at our hospital from 2008 to 2018. The pathologies for surgery were disc degeneration disease, spondylolisthesis, and segmental instability in the lumbar region and none of the selected patients had a history of previous spinal surgery. Based on a long-term follow-up study, the patients were divided into two groups [i.e., 1) ASD group (N = 15) and Non-ASD group (N = 15)]. The patients in the ASD and non-ASD groups were classified based on clinical and radiological indices after long-term follow-up investigations (5.37 ± 3.18 years). This study was approved by Chang Gung Memorial Hospital’s research ethics committee (approval No. 201702031B0) and signed informed consent was acquired from all participants prior to their enrolment in the relevant clinical protocol.
Methodology of personalized pre-operative finite element (FE) modeling of the lumbosacral spine form from lateral and anterior posterior (AP) X-ray radiographs. (A) Parameters extraction, (B) Development of the vertebrae, (C) Development of the FE model.
A non-linear poro-hyperelastic FE model of the lumbosacral spine (L1-S1) was developed for each patient based on their extracted geometrical values from pre-operative (pre-op) X-ray images (Figure 1). Each FE model consists of bony parts (i.e., posterior bony elements and vertebral bodies including cancellous and cortical bones) and soft tissues (i.e., five IVDs and ten cartilaginous endplates, seven ligaments, and five pairs of cartilaginous facet joints). The IVDs were characterized in FE modeling by a reinforced complex material consisting of the annulus fibrosus (AF) ground matrix reinforced with AF collagen fibers, and nucleus pulposus (NP). The drained solid phase of bony parts was considered isotropic elastic. However, the drained solid phase of the AF matrix and NP region were replicated based on the non-linear Mooney–Rivlin hyperelastic elastic-plastic hardening theory based on relevant studies in the literature (Schmidt et al., 2007; El-Rich et al., 2009). In addition, the theory of poroelastic (Argoubi and Shirazi-Adl, 1996; Ferguson et al., 2004) was reflected in the time-dependent response of the bony parts, cartilaginous endplates, and IVDs in the FE model. For this purpose, the values of permeability were reflected by variables based on the calculated void ratio in simulations based on the following equation (Argoubi and Shirazi-Adl, 1996; Ferguson et al., 2004),k=k0e1+e0e01+e2expM1+e1+e0−1Where k0 is the input initial permeability and e is defined as follows,e=∅f1−∅fWhere Øf is the porosity of the tissue which varies with material deformation during FE calculations. The composite structure of AF was mimicked by embedding six concentric reinforced fiber lamellae with an orientation of ±35° within a distance of 1 mm in the AF ground substance (Naserkhaki et al., 2016a). A constant boundary pore pressure equal to 0.25 MPa was imposed on all external surfaces of the IVDs as an additional constraint to include the swelling sensation in IVDs (Schmidt et al., 2010; Galbusera et al., 2011a). Ligaments were represented in the FE model using non-linear truss elements which were connected to bony parts and their length could be updated based on the patient’s geometrical input values. The mechanical properties of the ligaments were non-linear elastic based on available data in the literature and they could only be activated in tension (Shirazi-Adl et al., 1986a; Pintar et al., 1992). The mechanical properties of different components in this lumbosacral spine FE model are presented in Table 1.
Mechanical properties of the geometrical personalized poroelastic finite element model.
Component
Mechanical property behavior
Values of the mechanical properties
References
Cortical bone
Linear poroelastic
E = 12,000 MPa, ν = 0.3, k0 = 1 × 10−20 (m4/Ns), e = 0.02
Argoubi and Shirazi-Adl (1996),Goto et al. (2003),Ferguson et al. (2004),Schmidt et al. (2010),Galbusera et al. (2011b),Park et al. (2013)
Cancellous bone
Linear poroelastic
E = 200 MPa, ν = 0.25, k0 = 1 × 10−13 (m4/Ns), e = 0.4
Argoubi and Shirazi-Adl (1996),Ferguson et al. (2004),Schmidt et al. (2007),Schmidt et al. (2010),Galbusera et al. (2011b),Shih et al. (2013)
Endplate
Linear poroelastic
E = 5 MPa, ν = 0.1, k0 = 7.5 × 10−15 (m4/Ns), e = 4
Argoubi and Shirazi-Adl (1996),Goto et al. (2003),Ferguson et al. (2004),Schmidt et al. (2007),Schmidt et al. (2010),Galbusera et al. (2011b)
The vertebrae, cartilaginous endplates, and IVDs in each level were attached together in their anatomical positions using surface-to-surface tie contact conditions which provide equal translational and rotational motions at connected nodes. The facet joint surfaces were approximated in the model by a plane in which its orientation was defined by two card angles (Van Schaik et al., 1985; Panjabi et al., 1993). Based on the previous data in the literature (Panjabi et al., 1993) and our sensitivity analyses, the card angle about the x-axis was considered constant (80 degrees) and the one about the y-axis was chosen as the variable parameter which is extracted from the AP image. A surface-to-surface contact algorithm for both normal and tangential directions was considered to represent the articulation of the facet joints. A soft frictionless contact within an initial gap length of 0.5 mm to mimic the articulation of the facet joints in the FE model was considered (Naserkhaki et al., 2016a; Naserkhaki et al., 2016b; Naserkhaki and El-Rich, 2017). The transmitted force through contacting surfaces was mimicked using an exponential pressure-overclosure from zero at the initial gap until the contact pressure reached 120 MPa (Schmidt et al., 2010; Naserkhaki et al., 2016a; Naserkhaki et al., 2016b). To optimize the spine stability under compression while it lacks muscles, the weight of the upper body was applied as a compressive load using the follower load technique in which the line of action followed the spine curvature and passed through the vertebral bodies’ centroids (Patwardhan et al., 1999; Shirazi-Adl and Parnianpour, 2000; Dreischarf et al., 2014; Naserkhaki et al., 2016a; Naserkhaki and El-Rich, 2017). This follower load was applied to the models using pre-compressed unidirectional springs inserted between the centroids of two adjacent vertebral bodies. The rotational moments in different directions (i.e., flexion, extension, left/right lateral bending, and left/right axial rotation) were applied to the superior surface of the lumbosacral spine (L1) and Dirichlet boundary settings were applied at the sacral region to inhibit any displacement/rotation in all degrees of freedom. The model verification was approved based on mesh sensitivity analyses, and the finalized FE models contain 186,325 elements.
The overall validity of this FE modeling was previously confirmed for both static (Nikkhoo et al., 2020) and cyclic loading (Nikkhoo et al., 2015; Khalaf and Nikkhoo, 2021) which was well-aligned with both numerical and experimental studies (Rohlmann et al., 2009; Dreischarf et al., 2011; Dreischarf et al., 2012; Dreischarf et al., 2014). To evaluate the usability of the developed FE modeling technique in clinical applications, both pre- and post-operation functional X-ray images (patient in neutral, flexion, and extension positions) were employed to validate the predicted intersegmental range of motion (ROM). The functional X-ray images were used to measure the total lordosis angles (L1-S1) in neutral, flexion, and extension positions, while the ROM for each patient during flexion and extension was measured. The calculated subject-specific rotation of the L1 vertebra was then applied to the FE model, and the predicted rotations of each vertebra (i.e., intersegmental ROM) were compared with the measured ones from the images. To minimize the simplification errors regarding the consistency of the boundary conditions of the model with the in-vivo lumbar spine, a rotational control technique was utilized in the validation phase. To compare the achieved results, the percentage of the root mean square error (RMSE) was calculated as follows,RMSE=∑i=1i=NROMX−RayImagei−ROMFEMiROMX−RayImagei2Nwhere, “N” was considered equal to 5, corresponding to 5 spinal levels (L1-L2, to L5-S1) for the pre-op FE model.
2.2 Geometrical personalized FE modeling of the post-operative lumbosacral spine
To evaluate if the developed geometrical personalized FE models can predict any differences between biomechanical responses of the non-ASD group versus the ASD patient group, the post-operative (post-op) models of each patient were developed based on the X-ray images which were obtained 3 months post-surgery for both groups. To mimic the fusion surgery, a wide laminectomy was simulated by removing the relevant bony parts, IVD, posterior longitudinal ligament (PLL), and ligamentum flavum (LF). A posterior bilateral pedicle screw fixation system (including four pedicle screws, and two Titanium rods) was implanted and two posterior lumbar interbody fusion (PLIF) cages were inserted in the fusion level. The geometrical parameters of the fusion level (i.e., pedicle screw size, interbody cage height, and lordosis angle) were carefully adapted based on the post-op X-ray images. The material properties of the screws, rods, and interbody PLIF cages were considered isotropic elastic from literature (Zhang et al., 2018c) (Table 1). The tie contact boundary condition was set to constrain equal displacement during calculations (i.e., the same rotational and translational movement) for connected surfaces between the vertebrae, screws, and rods for mimicking the permanent fusion. In addition, the rod surfaces were directly connected to the screws using a tie contact technique to mimic fastening the screws using nuts in the fixation system.
Similar to the pre-op FE models, the accuracy of the predictions from post-op FE models were compared with those measured from functional post-op X-ray images during the movement in the sagittal plane. Incidentally, we removed the fused level from RMSE calculations for the post-op models, as this level was fixed during rotation, hence “N” was considered equal to four in Eq. 3. The post-op models for both the ASD and non-ASD groups were simulated under combined loading (i.e., compressive load and 10 Nm rotation in different planes) (Dreischarf et al., 2014) under the same aforementioned boundary conditions. To evaluate the time-dependent responses of the models subjected to cyclic loading, a daily cyclic loading scenario was applied to the FE models [i.e., 16 h of cyclic compressive loading of 500–1,000 N (40 and 20 min, respectively) after an 8 h pre-conditioning resting phase of 200 N (Nikkhoo et al., 2021)]. Different rotational movements (i.e., flexion, extension, right and left lateral bending, right and left axial rotation) were superimposed using a 10 N m moment before and after cyclic loading (i.e., points 1 and 2 in Figure 2) to model the rotational motions in the morning and the evening. The rotational moments were linearly applied and removed after 10s and only one motion was evaluated in each diurnal loading simulation. The biomechanical responses of the lumbosacral FE spine models in both groups were analyzed and compared before and after daily loading.
Loading scenario of the daily cyclic loading (Flexion, extension, lateral bending, and axial rotation moments of 10 N m were applied at points 1 and 2).
2.3 Statistical analyses for comparison of the results for ASD versus Non-ASD patients
The results of FE models (including the ROM, IVD height loss, intradiscal fluid loss, stress in AF matrix, and AF collagen fiber strain) in ASD and Non-ASD patients were compared. We used the non-parametric Friedman comparative test to specify the differences of the predicted results for ASD versus Non-ASD patients. The p-values less than 0.05 were reflected as significant statistical differences in this study.
3 Results
The geometrical personalized pre-op and post-op lumbosacral models were successfully developed for all thirty patients in both ASD and Non-ASD groups using our modeling updating algorithm and the relevant mesh sensitivity analyses confirmed the verification of their predictions. The percentage of calculated RMSE for pre-op FE models were 14.85% and 16.19% in flexion and 17.93% and 17.27% in extension for non-ASD and ASD patients’ models, respectively (Table 2). Correspondingly, these values for post-op models were 23.33% and 20.87% in flexion and 26.23% and 27.64% in extension for non-ASD and ASD patients, respectively (Table 3). Figure 3 schematically presents a sample of FE results overlapped with X-ray images for a patient from the non-ASD group.
Calculated percentage of root mean square errors (RMSE) for evaluating the accuracy and validity of the personalized pre-operative (pre-op) FE models for non-ASD patients (N = 15) and ASD patients (N = 15).
FE model
Percentage of the RMSE for measured and calculated ROM
Non-ASD group
ASD group
Flexion (%)
Extension (%)
Flexion (%)
Extension (%)
Patient No. 1
13.58
16.32
18.04
15.36
Patient No. 2
11.29
13.39
8.51
14.59
Patient No. 3
21.08
26.49
25.78
28.91
Patient No. 4
8.49
7.29
11.30
12.24
Patient No. 5
24.03
25.30
15.73
23.68
Patient No. 6
17.89
21.37
18.08
25.63
Patient No. 7
11.03
13.58
15.61
19.34
Patient No. 8
14.31
13.78
9.08
8.72
Patient No. 9
24.78
27.50
34.23
27.59
Patient No. 10
6.25
11.33
20.54
26.09
Patient No. 11
16.02
24.11
7.83
16.49
Patient No. 12
12.32
15.98
16.29
13.38
Patient No. 13
7.43
12.63
13.19
18.06
Patient No. 14
15.16
21.22
10.24
14.41
Patient No. 15
19.01
18.59
18.35
24.58
Average
14.85%
17.93%
16.19%
19.27%
Calculated percentage of root mean square errors (RMSE) for evaluating the accuracy and validity of the personalized post-operative (post-op) FE modeling technique.
FE model
Percentage of the RMSE for measured and calculated ROM
Non-ASD group
ASD group
Flexion (%)
Extension (%)
Flexion (%)
Extension (%)
Model No. 1
16.39
15.78
23.47
30.32
Model No. 2
17.84
22.07
12.49
25.28
Model No. 3
15.43
20.57
34.68
44.30
Model No. 4
12.48
17.42
19.66
23.79
Model No. 5
26.70
35.23
22.08
29.48
Model No. 6
18.48
19.54
15.41
24.54
Model No. 7
16.85
23.93
18.83
24.69
Model No. 8
22.73
28.02
8.40
21.18
Model No. 9
37.01
41.04
38.24
46.32
Model No. 10
9.20
19.07
28.74
33.75
Model No. 11
26.39
38.56
12.55
13.43
Model No. 12
19.69
24.23
21.35
30.61
Model No. 13
13.52
20.63
18.45
18.29
Model No. 14
38.67
36.38
17.83
21.04
Model No. 15
28.57
31.03
24.00
27.59
Average
23.33%
26.23%
20.87%
27.64%
A schematic sample of comparison between the finite element results and the functional X-ray images in the neutral position, flexion, and extension for (A) pre-operative and (B) post-operative models.
During the static loading scenario, the average ROM at the upper and lower adjacent levels increased post-surgery for both the non-ASD and ASD groups (Figure 4). The differences in average ROMs between pre-op and post-op results were significant for flexion and extension movements for both groups. However, there was no significant difference when we compared the values between the two groups. The significantly higher ROM in lateral bending was only observed in the ASD group for the upper adjacent level (Figure 4). The calculated values for IDP were quite similar in that significant differences were detected for flexion, extension, and lateral bending in both upper and lower adjacent levels (Figure 5). The FJF values showed a significant increase for extension and lateral bending as well (Figure 6).
Intersegmental motion patterns for non-ASD and ASD group FE models in the (A) upper and (B) lower adjacent level for different directions.
Calculated intradiscal pressure for non-ASD and ASD group FE models in the (A) upper and (B) lower adjacent level for different directions.
Calculated facet joint forces for non-ASD and ASD group FE models in the (A) upper and (B) lower adjacent level for different directions.
After daily cyclic loading (i.e., 16 h of cyclic loading to simulate the daily activities), the adjacent disc height averagely decreased by 5.56% and 7.17% in the pre-op models of non-ASD and ASD patients, respectively. These values were 11.03% and 16.03% for the post-op models in the non-ASD and ASD groups, respectively. Significant differences were observed between pre-op and post-op models in both groups (Figure 7). In addition, the differences between the post-op models in the non-ASD and ASD groups were significant for both upper and lower adjacent levels (p values equal to 0.016 and 0.021, respectively) (Figure 7). The fluid loss values were on average 13.44% and 15.29% for the pre-op models and 20.38% and 25.36% for the post-op models in the two groups after daily cyclic loading, respectively. The comparative difference between the post-op models in the non-ASD and ASD groups was significant (p-value = 0.039) for the upper adjacent level (Figure 7).
Disc height loss and fluid loss for post-operative non-ASD and ASD group FE models for the (A) upper and (B) lower adjacent levels.
The axial stress in the AF matrix significantly increased after fusion in sagittal plane movement (i.e., flexion and extension) and the differences between the non-ASD and ASD groups were significant in this plane only for the post-op models (Figure 8). However, fusion surgery did not significantly alter the AF axial stress in either lateral bending or axial rotation (Figure 8). Similar trends were calculated for collagen fiber strains. However, the differences between the non-ASD and ASD groups were significant for flexion and extension for both the pre-op and post-op results (Figure 9).
Increased axial stress in AF matrix for post-operative non-ASD and ASD group FE models for the (A) upper and (B) lower adjacent levels.
Increased collagen fiber strain in AF matrix for post-operative non-ASD and ASD group FE models for the (A) upper and (B) lower adjacent levels.
4 Discussion
The rigid instrumented posterolateral fusion (PLF) and posterior lumbar interbody fusion (PLIF) techniques are commonly considered the gold standard surgical treatment for different degenerative lumbar pathologies. Understanding lumbar spine biomechanics in post-surgical maneuvers could be beneficial for pre-planning the treatments based on primary rough estimations using personalized FE calculations. Hence, this study utilized a validated geometrically patient-specific FE modeling technique which can be used as a simple and cost-effective tool to evaluate the alteration of biomechanical response in adjacent segments post-fusion. Thirty patients were categorized for evaluation in this study in two distinct groups [i.e., (Serhan et al., 2011) non-ASD (N = 15) and (Kim and Choi, 2018) ASD patients (N = 15)] based on long-term clinical follow-up investigations. Although the current personalized FE modeling technique was previously validated against experimental in-vitro and numerical studies, the accuracy of each FE model for these 30 patients was evaluated based on the pre-op and post-op functional X-ray images in the sagittal plane. The achieved errors were averagely below 20% and 25% for pre-op and post-op models which confirms the applicability of this predictive algorithm. The observed error can be tolerated as this technique may be used as a rough estimation for clinicians. Hence, we developed and simulated pre-op and post-op FE models for each patient (in total, 60 FE models were simulated in this study) and evaluated the variations in biomechanical responses of adjacent segments post-surgery. We used the post-op X-ray images after 3 months in which no signs of ASD were observed to investigate if the results of FE simulations could show any significant difference between the non-ASD and ASD groups. Repeating the simulations based on acceptable numbers of patients (which was 15 patients in each group) included the influences of the patient’s anatomical parameters (such as disc height, shape and geometry of vertebra, lumbar lordotic angle before and after surgery, interbody cage height, etc.) to evaluate the aforementioned hypothesis of this study.
ASD may possibly be initiated and developed in 2–5 years based on induced modifications post-surgery, however, it is not possible to simulate this long-term phenomenon using numerical modeling. Meanwhile, most of the available FE studies which investigated the outcomes of different fixation surgical techniques only used static loading and evaluated the model response under simplified loading conditions (Jin et al., 2012; Jahng et al., 2013; Guo et al., 2019; Guo and Yin, 2019). Considering the time-dependent behavior of both vertebra and IVD in lumbar spine motion segments can enhance the accuracy of the biomechanical predictions. One of our contributions to the literature in this study was to utilize the non-linear poroelasticity theory in FE simulations to represent the interaction of interstitial water in a saturated solid matrix which makes the calculations more complicated and laborious. Hence, we provided more realistic simulations to evaluate the biomechanical responses of the lumbosacral spine under a daily cyclic loading regime (i.e., 8 h of resting time prior to 16 h of daily activity) using the aforementioned time-dependent constitutive equations. Using the poroelasticity theory, the interstitial water flow was calculated during daily cyclic loading to extract the variations in disc height and fluid loss. Although a daily cyclic loading (24 h) cannot represent a long-term loading condition, significant variations in output parameters (such as disc height, disc fluid loss, stress in the AF matrix, and strain in collagen fibers) may be characterized as indicators to reflect the initiation of abnormality in the biomechanical response which could be accumulated by repetitive loading.
The primary findings of this study demonstrated that the motion patterns in adjacent segments altered in post-op models. Although a slightly increased in ROM was calculated in most of the post-op FE models in different movement directions, the variation of ROM was significant in the sagittal plane (i.e., flexion and extension). However, no significant difference was observed between the calculated ROM of adjacent segments for the non-ASD and ASD groups. Posterior fusion results in increasing the spinal segmental rigidity which could possibly result in rotational compensation to the adjacent segments (Chiang et al., 2006; Sudo et al., 2006). Based on this hypothesis, the IDP and FJF values were significantly increased in adjacent levels in extension and lateral bending movements as well. Similar to ROM, this static loading condition did not demonstrate any significant difference between the non-ASD and ASD patients.
Simulating the biomechanical response of the pre-op and post-op lumbosacral spine FE models subjected to the daily cyclic loading provided an enhanced possibility to include the effect of damping characterization in the achieved results. For this purpose, disc height loss and fluid loss were chosen in the current comparative investigation as two critical indicators of the initiation of IVD denaturation which may lead to mild/severe IVD degeneration in a couple of years. Applying repetitive compressive loading to pre-op models results in an approximately uniform disc height reduction which is consistent with previous studies in the literature (Tyrrell et al., 1985; Lai et al., 2008). The induced modification regarding the interbody fusion in post-op models alters the load sharing and motion patterns by increasing the stiffness of the instrumented segment which leads to increased IVD height loss and fluid loss. Although no statistical differences were calculated for pre-op FE results between the two groups, the post-op FE results demonstrated that both IVD height loss and fluid loss significantly increased for the patients with ASD. Hence, the results of this study highlighted the effect of geometrical parameters (which may refer to the anatomical conditions or the induced modifications regarding surgical techniques) on the time-dependent response of lumbar spine biomechanics. The variations in IVD height loss and fluid loss may perhaps indicate initiation of ASD as it is confirmed for more than half of patients who underwent rigid-fusion surgery in previous clinical investigations (Miyakoshi et al., 2000; Ishihara et al., 2001). The alteration of fluid-solid interaction during the cyclic loading possibly changes the stress/strain distribution in the AF matrix. Increasing the disc height loss and fluid loss decreases the contribution of the fluid resistance in IVD which leads to an increase of stress in the AF matrix and collagen fibers. This is the reason that no significant differences were calculated for IDP between the two groups in static loading but the stress/strain increased after 16 h of cyclic loading. Analogous biomechanical response patterns were calculated for the stress and fiber strain in the disc AF matrix region post-fusion. Increased axial stress was significant in both groups when comparing the pre-op and post-op results for the sagittal plane. However, the ASD group showed significantly higher stress for both flexion and extension movement. An increased amount of stress in the AF matrix refers to decreasing the effect of damping resistive interstitial water in total IVD bulk stiffness (Schmidt et al., 2010; Galbusera et al., 2011a). Hence, it can be related to a higher rate of fluid loss rate after repetitive daily loading. On the other hand, rigid fusion surgery can alter the load sharing by removing the load path from the anterior region (vertebral bodies and IVDs) and shifting it to the posterior fixation system which significantly affects the adjacent IVD loading conditions. In addition, the increased fiber strain in the AF region was observed in upper adjacent levels between the two groups for both the pre-op and post-op models which clearly highlights that this geometrically patient-specific FE modeling can approximately differentiate the non-ASD patients from the patients with ASD.
Some limitations of this work should be acknowledged. We used a geometrically personalized modeling algorithm which means that this methodology is only sensitive to anatomical conditions and the same mechanical properties were assumed for the different components of the lumbosacral spine FE models for both groups. As there is no feasible algorithm to extract the mechanical properties of different tissues from X-ray images and it is very difficult to develop optimization methodology for this purpose, the observed variations in motion patterns and load sharing for different patients are only based on different geometrical conditions in this study. In addition, the main contribution of this study is to investigate if the anatomical parameters alter the biomechanical responses of adjacent segments following lumbar fusion surgery. Hence, this simplification might be tolerated to make this study feasible. Another limitation is neglecting the effect of active muscle force in this geometrically personalized osseoligamentous FE modeling technique. Although the follower load approach (Patwardhan et al., 1999; Shirazi-Adl and Parnianpour, 2000; Dreischarf et al., 2014) was considered in this study to mimic the passive response of muscles and upper body weight, there is a lack of muscle force reaction which may alter the biomechanical response specially for possible post-surgical muscle damages. To augment the personalized muscle forces to the model, new sets of clinical measurements before and after surgery are needed which can be planned in our future works.
5 Conclusion
This study used a validated geometrically personalized FE modeling algorithm which has the potential to be used in clinics to analyze the biomechanical response of adjacent segments post posterior lumbar fusion surgery. Evaluating the biomechanical response of pre-op and post-op modeling in the non-ASD versus ASD groups showed that the geometrical differences among patients (such as vertebral geometries, IVD height, lumbar lordosis angle, interbody cage size) cause significant variations in the estimated mechanical response. In addition, significant differences between estimated disc height loss, fluid loss, stress, and strain in adjacent levels for aforementioned patient groups can predict the increased risk of pathological development of ASD, which is consistent with our long-term clinical observation in this study. In conclusion, the results of the current study highlighted the effect of geometrical parameters (which may refer to the anatomical conditions or the induced modifications regarding surgical techniques) on the time-dependent response of lumbar spine biomechanics. Henceforward, performing the simulations based on geometrical personalized FE models which include the patient’s anatomical parameters can be used as a surgical planning tool in clinical settings to minimize further long-term complications. It may provide clinicians with a valuable pre-planning tool to make informed pre-op and post-op decisions.
Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors upon request.
Ethics statement
The studies involving human participants were reviewed and approved by Chang Gung Memorial Hospital’s research ethics committee (approval No. 201702031B0). The patients/participants provided their written informed consent to participate in this study.
Author contributions
MN, W-CC, M-LL, C-JF, C-CN, H-YL, and C-HC conceived and designed this study. W-CC, M-LL, C-JF, and C-CN provided the clinical study. MN and C-HC conducted FE modeling and statistical analyses. MN, W-CC, H-YL, and C-HC prepared the manuscript. MN, W-CC, M-LL, C-JF, C-CN, H-YL, and C-HC commented, revised the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding
The authors acknowledge the funding supported by the National Science and Technology Council of Taiwan (111-2221-E-182-009-MY3), the Chang Gung Memorial Hospital Research Program (CMRPD1L0181, CMRPD1L0182, CMRPG3M0481), and the Healthy Aging Research Center, Chang Gung University, Taoyuan, Taiwan (EMRPD1M0411).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
ReferencesArgoubiM.Shirazi-AdlA. (1996). Poroelastic creep response analysis of a lumbar motion segment in compression. J. Biomechanics29 (10), 1331–1339. 10.1016/0021-9290(96)00035-8BeckmannA.NicoliniL. F.GrevensteinD.BackesH.OikonomidisS.SobottkeR. (2020). Biomechanical in vitro test of a novel dynamic spinal stabilization system incorporating polycarbonate urethane material under physiological conditions. J. Biomechanical Eng.142 (1), 011005. 10.1115/1.4044242ChiangM. F.ZhongZ. C.ChenC. S.ChengC. K.ShihS. L. (2006). Biomechanical comparison of instrumented posterior lumbar interbody fusion with one or two cages by finite element analysis. Spine (Phila Pa 1976)31 (19), E682–E689. 10.1097/01.brs.0000232714.72699.8eDreischarfM.RohlmannA.BergmannG.ZanderT. (2012). Optimised in vitro applicable loads for the simulation of lateral bending in the lumbar spine. Med. Eng. Phys.34 (6), 777–780. 10.1016/j.medengphy.2012.04.002DreischarfM.RohlmannA.BergmannG.ZanderT. (2011). Optimised loads for the simulation of axial rotation in the lumbar spine. J. Biomech.44 (12), 2323–2327. 10.1016/j.jbiomech.2011.05.040DreischarfM.ZanderT.Shirazi-AdlA.PuttlitzC. M.AdamC. J.ChenC. S. (2014). Comparison of eight published static finite element models of the intact lumbar spine: Predictive power of models improves when combined together. J. Biomech.47 (8), 1757–1766. 10.1016/j.jbiomech.2014.04.002EbrahimkhaniM.ArjmandN.Shirazi-AdlA. (2022). Adjacent segments biomechanics following lumbar fusion surgery: A musculoskeletal finite element model study. Eur. Spine J.31 (7), 1630–1639. 10.1007/s00586-022-07262-3El-RichM.ArnouxP. J.WagnacE.BrunetC.AubinC. E. (2009). Finite element investigation of the loading rate effect on the spinal load-sharing changes under impact conditions. J. Biomech.42 (9), 1252–1262. 10.1016/j.jbiomech.2009.03.036ErbulutD. U.ZafarparandehI.OzerA. F.GoelV. K. (2013). Biomechanics of posterior dynamic stabilization systems. Adv. Orthop.2013, 1–6. 10.1155/2013/451956FergusonS. J.ItoK.NolteL-P. (2004). Fluid flow and convective transport of solutes within the intervertebral disc. J. Biomechanics37 (2), 213–221. 10.1016/s0021-9290(03)00250-1GalbuseraF.SchmidtH.Neidlinger-WilkeC.WilkeH-J. (2011). The effect of degenerative morphological changes of the intervertebral disc on the lumbar spine biomechanics: A poroelastic finite element investigation. Comput. Methods Biomechanics Biomed. Eng.14 (8), 729–739. 10.1080/10255842.2010.493522GalbuseraF.SchmidtH.NoaillyJ.MalandrinoA.LacroixD.WilkeH-J. (2011). Comparison of four methods to simulate swelling in poroelastic finite element models of intervertebral discs. J. Mech. Behav. Biomed. Mater.4 (7), 1234–1241. 10.1016/j.jmbbm.2011.04.008GotoK.TajimaN.ChosaE.TotoribeK.KuboS.KurokiH. (2003). Effects of lumbar spinal fusion on the other lumbar intervertebral levels (three-dimensional finite element analysis). J. Orthop. Sci.8 (4), 577–584. 10.1007/s00776-003-0675-1GuoL. X.YinJ. Y. (2019). Finite element analysis and design of an interspinous device using topology optimization. Med. Biol. Eng. Comput.57 (1), 89–98. 10.1007/s11517-018-1838-8GuoT-M.LuJ.XingY-L.LiuG-X.ZhuH-Y.YangL. (2019). A 3-dimensional finite element analysis of adjacent segment disk degeneration induced by transforaminal lumbar interbody fusion after pedicle screw fixation. World Neurosurg.124, e51–e57. 10.1016/j.wneu.2018.11.195HelgesonM. D.BevevinoA. J.HilibrandA. S. (2013). Update on the evidence for adjacent segment degeneration and disease. Spine J.13 (3), 342–351. 10.1016/j.spinee.2012.12.009HsiehC-T.ChangC-J.SuI. C.LinL-Y. (2016). Clinical experiences of dynamic stabilizers: Dynesys and Dynesys top loading system for lumbar spine degenerative disease. Kaohsiung J. Med. Sci.32 (4), 207–215. 10.1016/j.kjms.2016.03.007IshiharaH.OsadaR.KanamoriM.KawaguchiY.OhmoriK.KimuraT. (2001). Minimum 10-year follow-up study of anterior lumbar interbody fusion for isthmic spondylolisthesis. J. Spinal Disord.14 (2), 91–99. 10.1097/00002517-200104000-00001JahngT-A.KimY. E.MoonK. Y. (2013). Comparison of the biomechanical effect of pedicle-based dynamic stabilization: A study using finite element analysis. Spine J.13 (1), 85–94. 10.1016/j.spinee.2012.11.014JinY. J.KimY. E.SeoJ. H.ChoiH. W.JahngT-A. (2012). Effects of rod stiffness and fusion mass on the adjacent segments after floating mono-segmental fusion: A study using finite element analysis. Eur. Spine J.22 (5), 1066–1077. 10.1007/s00586-012-2611-6KhalafK.NikkhooM. (2021). Comparative biomechanical analysis of rigid vs. flexible fixation devices for the lumbar spine: A geometrically patient-specific poroelastic finite element study. Comput. Methods Programs Biomed.212, 106481. 10.1016/j.cmpb.2021.106481KimJ. E.ChoiD. J. (2018). Biportal endoscopic transforaminal lumbar interbody fusion with arthroscopy. Clin. Orthop. Surg.10 (2), 248–252. 10.4055/cios.2018.10.2.248KimJ. Y.RyuD. S.PaikH. K.AhnS. S.KangM. S.KimK. H. (2016). Paraspinal muscle, facet joint, and disc problems: Risk factors for adjacent segment degeneration after lumbar fusion. Spine J.16 (7), 867–875. 10.1016/j.spinee.2016.03.010KimY. J.LeeS. G.ParkC. W.SonS.KimW. K. (2012). Long-term follow-up (minimum 5 Years) study of single-level posterior dynamic stabilization in lumbar degenerative disease; 'interspinous U' & 'DIAM. Korean J. Spine9 (2), 102. 10.14245/kjs.2012.9.2.102LaiA.ChowD. H.SiuS. W.LeungS. S.LauE. F.TangF. H. (2008). Effects of static compression with different loading magnitudes and durations on the intervertebral disc: An in vivo rat-tail study. Spine (Phila Pa 1976)33 (25), 2721–2727. 10.1097/brs.0b013e318180e688LavilleA.LaporteS.SkalliW. (2009). Parametric and subject-specific finite element modelling of the lower cervical spine. Influence of geometrical parameters on the motion patterns. J. Biomech.42 (10), 1409–1415. 10.1016/j.jbiomech.2009.04.007LiangJ.DongY.ZhaoH. (2014). Risk factors for predicting symptomatic adjacent segment degeneration requiring surgery in patients after posterior lumbar fusion. J. Orthop. Surg. Res.9, 97. 10.1186/s13018-014-0097-0MiyakoshiN.AbeE.ShimadaY.OkuyamaK.SuzukiT.SatoK. (2000). Outcome of one-level posterior lumbar interbody fusion for spondylolisthesis and postoperative intervertebral disc degeneration adjacent to the fusion. Spine25 (14), 1837–1842. 10.1097/00007632-200007150-00016NaserkhakiS.El-RichM. (2017). Sensitivity of lumbar spine response to follower load and flexion moment: Finite element study. Comput. Methods Biomech. Biomed. Engin20 (5), 550–557. 10.1080/10255842.2016.1257707NaserkhakiS.JaremkoJ. L.AdeebS.El-RichM. (2016). On the load-sharing along the ligamentous lumbosacral spine in flexed and extended postures: Finite element study. J. Biomech.49 (6), 974–982. 10.1016/j.jbiomech.2015.09.050NaserkhakiS.JaremkoJ. L.El-RichM. (2016). Effects of inter-individual lumbar spine geometry variation on load-sharing: Geometrically personalized Finite Element study. J. Biomech.49 (13), 2909–2917. 10.1016/j.jbiomech.2016.06.032NikkhooM.KhalafK.KuoY-W.HsuY-C.HaghpanahiM.ParnianpourM. (2015). Effect of degeneration on fluid-solid interaction within intervertebral disk under cyclic loading - a meta-model analysis of finite element simulations. Front. Bioeng. Biotechnol.3, 4. 10.3389/fbioe.2015.00004NikkhooM.KhozZ.ChengC-H.NiuC-C.El-RichM.KhalafK. (2020). Development of a novel geometrically-parametric patient-specific finite element model to investigate the effects of the lumbar lordosis angle on fusion surgery. J. Biomechanics102, 109722. 10.1016/j.jbiomech.2020.109722NikkhooM.LuM. L.ChenW. C.FuC. J.NiuC. C.LinY. H. (2021). Biomechanical investigation between rigid and semirigid posterolateral fixation during daily activities: Geometrically parametric poroelastic finite element analyses. Front. Bioeng. Biotechnol.9, 646079. 10.3389/fbioe.2021.646079PanjabiM. M.OxlandT.TakataK.GoelV.DuranceauJ.KragM. (1993). Articular facets of the human spine. Quantitative three-dimensional anatomy. Spine (Phila Pa 1976)18 (10), 1298–1310. 10.1097/00007632-199308000-00009ParkW. M.KimK.KimY. H. (2013). Effects of degenerated intervertebral discs on intersegmental rotations, intradiscal pressures, and facet joint forces of the whole lumbar spine. Comput. Biol. Med.43 (9), 1234–1240. 10.1016/j.compbiomed.2013.06.011PatwardhanA. G.HaveyR. M.MeadeK. P.LeeB.DunlapB. (1999). A follower load increases the load-carrying capacity of the lumbar spine in compression. Spine24 (10), 1003–1009. 10.1097/00007632-199905150-00014PhillipsF. M.VoronovL. I.GaitanisI. N.CarandangG.HaveyR. M.PatwardhanA. G. (2006). Biomechanics of posterior dynamic stabilizing device (DIAM) after facetectomy and discectomy. Spine J.6 (6), 714–722. 10.1016/j.spinee.2006.02.003PintarF. A.YoganandanN.MyersT.ElhagediabA.SancesA. (1992). Biomechanical properties of human lumbar spine ligaments. J. Biomechanics25 (11), 1351–1356. 10.1016/0021-9290(92)90290-hRohlmannA.ZanderT.RaoM.BergmannG. (2009). Realistic loading conditions for upper body bending. J. Biomech.42 (7), 884–890. 10.1016/j.jbiomech.2009.01.017SchmidtH.HeuerF.SimonU.KettlerA.RohlmannA.ClaesL. (2006). Application of a new calibration method for a three-dimensional finite element model of a human lumbar annulus fibrosus. Clin. Biomech. (Bristol, Avon)21 (4), 337–344. 10.1016/j.clinbiomech.2005.12.001SchmidtH.KettlerA.HeuerF.SimonU.ClaesL.WilkeH. J. (2007). Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading. Spine (Phila Pa 1976)32 (7), 748–755. 10.1097/01.brs.0000259059.90430.c2SchmidtH.Shirazi-AdlA.GalbuseraF.WilkeH. J. (2010). Response analysis of the lumbar spine during regular daily activities--a finite element analysis. J. Biomech.43 (10), 1849–1856. 10.1016/j.jbiomech.2010.03.035SerhanH.MhatreD.DefossezH.BonoC. M. (2011). Motion-preserving technologies for degenerative lumbar spine: The past, present, and future horizons. SAS J.5 (3), 75–89. 10.1016/j.esas.2011.05.001ShihS. L.LiuC. L.HuangL. Y.HuangC. H.ChenC. S. (2013). Effects of cord pretension and stiffness of the Dynesys system spacer on the biomechanics of spinal decompression-a finite element study. BMC Musculoskelet. Disord.14, 191. 10.1186/1471-2474-14-191Shirazi AdlA.AhmedA. M.ShrivastavaS. C. (1986a). Mechanical response of a lumbar motion segment in axial torque alone and combined with compression. Spine11 (9), 914–927. 10.1097/00007632-198611000-00012Shirazi-AdlA.AhmedA. M.ShrivastavaS. C. (1986b). A finite element study of a lumbar motion segment subjected to pure sagittal plane moments. J. Biomechanics19 (4), 331–350. 10.1016/0021-9290(86)90009-6Shirazi-AdlA.ParnianpourM. (2000). Load-bearing and stress analysis of the human spine under a novel wrapping compression loading. Clin. Biomech.15 (10), 718–725. 10.1016/s0268-0033(00)00045-0SudoH.OdaI.AbumiK.ItoM.KotaniY.MinamiA. (2006). Biomechanical study on the effect of five different lumbar reconstruction techniques on adjacent-level intradiscal pressure and lamina strain. J. Neurosurg. Spine5 (2), 150–155. 10.3171/spi.2006.5.2.150TyrrellA. R.ReillyT.TroupJ. D. (1985). Circadian variation in stature and the effects of spinal loading. Spine (Phila Pa 1976)10 (2), 161–164. 10.1097/00007632-198503000-00011Van SchaikJ. P.VerbiestH.Van SchaikF. D. (1985). The orientation of laminae and facet joints in the lower lumbar spine. Spine (Phila Pa 1976)10 (1), 59–63. 10.1097/00007632-198501000-00009VergroesenP. P.KingmaI.EmanuelK. S.HoogendoornR. J.WeltingT. J.van RoyenB. J. (2015). Mechanics and biology in intervertebral disc degeneration: A vicious circle. Osteoarthr. Cartil.23 (7), 1057–1070. 10.1016/j.joca.2015.03.028ZhangY.ZhangZ-C.LiF.SunT-S.ShanJ-L.GuanK. (2018). Long-term outcome of dynesys dynamic stabilization for lumbar spinal stenosis. Chin. Med. J.131 (21), 2537–2543. 10.4103/0366-6999.244107ZhangZ.FogelG. R.LiaoZ.SunY.LiuW. (2018). Biomechanical analysis of lumbar interbody fusion cages with various lordotic angles: A finite element study. Comput. Methods Biomech. Biomed. Engin21 (3), 247–254. 10.1080/10255842.2018.1442443ZhangZ.LiH.FogelG. R.LiaoZ.LiY.LiuW. (2018). Biomechanical analysis of porous additive manufactured cages for lateral lumbar interbody fusion: A finite element analysis. World Neurosurg.111, e581–e591. 10.1016/j.wneu.2017.12.127ZhaoX.DuL.XieY.ZhaoJ. (2018). Effect of lumbar lordosis on the adjacent segment in transforaminal lumbar interbody fusion: A finite element analysis. World Neurosurg.114, e114–e120. 10.1016/j.wneu.2018.02.073