Biomechanical Analysis of Adjacent Segments after Spinal Fusion Surgery Using a Geometrically Parametric Patient-Specific Finite Element Model

Abstract

This study aimed to perform a mechanical analysis of adjacent segments after spinal fusion surgery using a geometrically parametric patient-specific finite element model to elucidate the mechanism of adjacent segment degeneration (ASD), thereby providing theoretical evidence for early disease prevention. Fourteen parameters based on patient-specific spinal geometry were extracted from a patient's preoperative computed tomography (CT) scan, and the relative positions of each spinal segment were determined using the image match method. A preoperative patient-specific model of the spine was established through the above method. The postoperative model after L4-L5 posterior lumbar interbody fusion (PLIF) surgery was constructed using the same method except that the lamina and intervertebral disc were removed, and a cage, 4 pedicle screws, and 2 connecting rods were inserted. Range of motion (ROM) and stress changes were determined by comparing the values of each anatomical structure between the preoperative and postoperative models. The overall ROM of the lumbar spine decreased after fusion, while the ROM, stress in the facet joints, and stress in the intervertebral disc of adjacent segments all increased. An analysis of the stress distribution in the annulus fibrosus, nucleus pulposus, and facet joints also showed that not only was the maximum stress in these tissues elevated, but the areas of moderate-to-high stress were also expanded. During torsion, the stress in the facet joints and annulus fibrosus of the proximal adjacent segment (L3-L4) increased to a larger extent than that in the distal adjacent segment (L5-S1). While fusion surgery causes an overall restriction of motion in the lumbar spine, it also causes more load sharing by the adjacent segments to compensate for the fused segment, thus increasing the risk of ASD. The proximal adjacent segment is more prone to degeneration than the distal adjacent segment after spinal fusion due to the significant increase in stress.

Introduction

Intervertebral spinal fusion surgery is the most commonly used surgical procedure for the treatment of degenerative diseases of the lumbar spine¹. An excellent outcome in the short-term period after surgery can be achieved for more than 90% of patients². However, the results of a long-term follow-up study revealed that some patients developed degeneration of segments adjacent to the fused segment³. Lumbar interbody fusion accelerates degenerative changes in adjacent segments, which is known as adjacent segment degeneration (ASD). According to the literature, the incidence of ASD diagnosed based on medical imaging examinations ranges from 36% to 84% five years after fusion surgery⁴, which could lead to symptoms such as radiating pain or intermittent claudication and possibly even the need for revision surgery. The mechanism of ASD remains unknown, but most researchers believe that biomechanical factors play an important role. Some have attributed ASD to increased range of motion (ROM) of the adjacent segments after surgery^5,6, some have attributed it to increased intradiscal pressure in the adjacent segments^7,8,9, and others have attributed it to increased stress in the facet joints of the adjacent segments¹⁰.

Among the various methods used to study spinal biomechanics, finite element (FE) modeling is widely used because it is noninvasive, inexpensive, and reproducible. Some researchers^11,12,13 have established a 3D FE model of the whole lumbar spine (L1-L5) with data extracted from preoperative computed tomography (CT) scans, which made it possible to explore various aspects of spinal biomechanics, ranging from the response of the spine to different loading conditions^14,15 to the effects of different pathologies¹⁶ and the effects of relevant treatment modalities and techniques¹⁷. Although the above modeling method could provide output regarding the patient-specific geometry of the spine with a complex interface and a wealth of information otherwise unattainable from in vivo experiments, its clinical use has remained limited due to the time-consuming nature of the process, rendering the method available only for models based on one or a few subjects¹⁴. To address this problem, Nikkhoo et al.¹⁸ established a simplified L1-S1 lumbosacral model in which the geometry of the spine is controlled by parameters extracted from patients' preoperative image data, allowing patient-specific models to be automatically generated or updated according to the input parameters. The FE model based on this modeling method has been proven to have good validity. However, there were significant differences in the intradiscal pressure, mean stresses in the facet joints, and mean stresses in the annulus fibrosus compared with the previous CT-based reconstructed model. Another simplified spinal model was applied in a study by Ghezelbash et al.¹⁹, but this model differed greatly from the real geometry of the lumbar spine due to the cylindrical shape of the vertebrae and lack of structure regarding the posterior elements.

Therefore, in this study, we developed a geometrically parametric patient-specific FE model to achieve a more efficient modeling and analysis process with good validity. Then, we performed a mechanical analysis of adjacent segments after fusion surgery to elucidate the mechanism and provide theoretical evidence for the early prevention of ASD.

Protocol

The protocol was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Institutional Review Board of the China-Japan Friendship Hospital.

1. Parametric modeling of lumbar spine geometry

1. Extract the initial data (DICOM 3.0 format with a pixel size of 0.33 mm and a layer spacing of 1 mm) for modeling from a CT scan data set of an adult healthy male with no history of trauma, deformity or tumor of the spine (height 180 cm, weight 68 kg).
2. Select 14 characteristic parameters to achieve spinal contour generation, considering the morphological features of the lumbar spine that are most concerning in clinical practice and the latest literature^18,20.
   1. Measure all 14 of these parameters directly on CT images using 3D image processing software, as shown in Figure 1A.
   2. In the Axial View window, employ the Ellipse Tool to precisely measure the parameters of the vertebral endplate.
   3. For each spinal segment and vertebra, initially, scrutinize the CT images in the axial direction from top to bottom. Identify the image showcasing the most complete and largest area of the vertebral boundary for subsequent data measurement.
   4. Utilize the Ellipse Tool within the measurement module to conform to the lower endplate of the vertebra, as depicted in Figure 1A.
   5. After conducting multiple measurements, compute the average value to ensure that the variance between the ellipse area and the vertebral area remains within 10%.
   6. Gauge the length of the long and short axes of the fitted ellipse, denoting them as parameters A1 and A2.
   7. Additionally, measure the cross-section at the narrowest part of the middle section and the upper endplate, designating them as parameters B1, B2, C1, and C2.
   8. In the Coronal View window, employ the Distance Measurement Tool to determine the vertebral height, represented by parameter H.
3. Utilize the Angle Measurement Tool to quantify the posterior tilt angles of the upper and lower facets in the sagittal view, identified as parameters α and β.
   1. In the CT Axial View window, use the Distance Measurement Tool to measure the vertical distance from the middle section of the vertebra to the lamina as the pedicle length control parameter L1.
   2. In a similar manner, use the Angle Measurement Tool and the Distance Measurement Tool to measure the parameters L2, γ, L3, and θ in the CT axial view window, respectively, as the control parameters for the transverse process and the spinous process.
   3. To control for interobserver and intraobserver variability, let two physicians with more than 5 years of training in spinal surgery measure each of the parameters 3 times to confirm the reliability of the data.
4. Apply modeling software to build the model using the Solid Release function according to the simplified model design scheme in Figure 1A for all spinal segments except the sacrum.
   1. Establish three reference planes and adjust the distance between the top and bottom planes to match the vertebral height. On each plane, sketch three concentric ellipses, aligning their dimensions with the CT data measurements.
   2. Utilize these sketched ellipses as constraint contours for the Solid Release function, resulting in the creation of a simplified vertebral model.
5. To recreate the arc surface contact at the facet joint, employ Cylindrical Surfaces to mimic the facet surfaces.
   1. Ensure that the upper facet takes the form of a concave 1/4 cylindrical arc surface while the lower facet takes the form of a convex 1/4 cylindrical arc surface. To mitigate stress concentration during facet alignment, appropriately round the edges of the upper and lower facets.
6. Generate an extruded entity between the facet and the vertebra to emulate the pedicle. Since the shape of the transverse process or the spinous process minimally impacts subsequent ligament element additions, use a regular parallelepiped to replicate the geometric contour of these processes.
   1. Round some of the corners for a refined representation. Modify the values of the 14 feature parameters in the modeling software to generate a patient-specific spinal geometry.
7. Measure only the parameters C1 and C2 of the upper endplate of the sacrum in the CT cross-sectional window and the upper facet inclination angle parameter α in the sagittal window with a similar measurement method as described in step 1.1, considering that the FE model calculation mainly focuses on the stress of the facet and the upper endplate.
   1. Generate a cone structure with a wide top and a narrow bottom as a simplified model of the sacrum, with a columnar-like structure extending from both sides to mimic the sacral wings and provide ligament attachment points in the modeling software. Use the three parameters above to control the S1 geometry. See Figure 1B for the simplified sacral model.
8. Apply the image-matching method to determine the relative positions of each spinal segment.
   1. Import all vertebral and sacral models into the assembly interface of the modeling software, where a middle sagittal view of the CT image is loaded as the reference background.
   2. Rotate, move, and scale each vertebral segment to match it with the corresponding part of the reference image (Figure 1B).
9. Extract contours of the adjacent vertebral endplate for solid release.
   1. Select the adjacent vertebral endplates and insert them into the sketch to obtain the intervertebral disc.
   2. Use the Convert Entity Reference command to extract the vertebral endplate contours as elliptical lines in the sketch and perform sketch release to generate a simplified disc matrix model.
   3. Create the nucleus pulposus in a similar manner as the disc matrix, in addition to reducing the ellipse sketch to 40% of the original area and moving it back slightly in the sketch.
   4. Additionally, enlarge the nucleus pulposus model by 10% to facilitate end plate segmentation when meshing. See Figure 1B for the final simplified disc model.

2. Construction of the posterior lumbar interbody fusion (PLIF) model with patient-specific geometry

1. Reload the simplified model in the modeling software. Select the L4-L5 intervertebral disc segment for fusion.
2. Remove the lamina and spinous process of the L4 vertebra manually based on the individualized parametric model of the lumbar spine. Remove the L4-L5 intervertebral disc.
3. Place a cage in the intervertebral space for bone fusion, and fill the remaining intervertebral space around the fusion cage with bone structure.
4. Apply the posterior lumbar intervertebral fusion (PLIF) method by inserting pedicle screws in the pedicle bilaterally.
   1. According to the literature²¹, use screws 5.5 mm in diameter and 45 mm in length, fixation rods 6 mm in diameter and 60 mm in length, and a graft cage 22 mm in length and 8 mm in width.
   2. Adjust the position of the pedicle screws such that the entry point is at approximately the center of the pedicle.
   3. Employ the Boolean operation method in modeling software, utilizing the combination command within the feature options.
   4. Set the operation type to subtract, utilizing the L4 and L5 vertebrae as the primary entities and the pedicle screws as the subtractive entities to accomplish modeling of the L4 and L5 pedicle screw trajectories.
   5. Employ the same procedure, setting the operation type to add, to consolidate the screw and fixation rod models into a unified whole. See Figure 1B for the constructed patient-specific PLIF model.

3. Establishment of parametric, patient-specific, preoperative, and postoperative FE models

1. Mesh generation
   1. Employ mesh software²² to mesh the preoperative and postoperative models after geometric processing. Import the stp model and utilize the 2D Meshing Auto Mesh module to set surface mesh sizes and element types. Generate the surface mesh of the models.
   2. Use the 3D Meshing Solid Map module to set the entity mesh element types and automatically generate the solid mesh. Employ quadrilateral elements 1 mm in size for surface meshing of pedicle screws and fixation rods and automatically generate a solid mesh with a mix of C3D4 and C3D8R elements.
   3. For the intervertebral fusion device, employ triangular elements 1 mm in size for surface meshing and use C3D4 tetrahedral elements for the solid mesh.
   4. Due to the irregular shape of the residual intervertebral disc model for the L4-L5 segment, use triangular elements 1.5 mm in size for surface meshing on end faces. Generate the solid mesh through extrusion, utilizing a mix of C3D8R and C3D4 elements.
   5. Mesh the remaining parts of the post-PLIF model with the same method as the preoperative model, resulting in the creation of 617,231 cells and 151,078 nodes in the post-PLIF model.
2. Material properties and interaction settings
   1. Import the meshed pre- and postoperative models into FE software for preprocessing.
      1. In the Material Manager panel, set the Material Behavior of the pedicle screws, fixation rods, and intervertebral fusion devices as isotropic linear elastic materials.
      2. In the Data tab, specify the Young Modulus and Poisson Ratio of the materials.
      3. For the screws and rods, use titanium alloy, and for the fusion device, use polyetheretherketone. Refer to Table 2 for the specific material property parameters of these two materials.
      4. As the mesh of the residual intervertebral disc at L4-L5 is not hexahedral and cannot be defined as a hyperelastic material, refer to the relevant literature²³ and set it as an isotropic linear elastic material in the Material Manager; specify its Young Modulus as 4 MPa and Poisson Ratio as 0.45.
   2. Navigate to the Interaction module, open the Constraint Manager, click on the Create button to open the Create Constraint window.
      1. Set the type as Binding. In the Model Display window, select the upper and lower vertebral end faces, as well as the tied nodes of the fusion device.
      2. After confirming, open the Edit Constraint window, set the Discretization Method to the analysis default, and specify not to exclude shell element thickness.
      3. Set the interaction settings according to the biomechanical conditions after ideal intervertebral fusion. Ignore possible slippage between the bone and the screws or cage.
      4. Set the contact relations between the screws and cancellous bone and between the cage and the end surfaces of the upper and lower vertebral bodies as binding.
      5. Set the contact interaction property between the joint contact surfaces as sliding friction controlled by the Penalty function, with a friction coefficient of 0.01 in the tangential direction and hard contact in the normal direction where separation after contact is allowed.
   3. Set the boundary conditions according to the rules of human lumbosacral motion, where all vertebral segments can move, while the sacrum mainly provides support and fixation.
      1. Access the Load module in FE software, open the Boundary Conditions Manager, and click on the Create button to open the Create Boundary Condition window.
      2. Set the Category as Mechanical and choose the type applicable to the selected analysis step as Symmetry/Anti-symmetry/Full fixity.
      3. Click Continue, and in the model display interface, select the surface nodes of the sacrum.
      4. Upon completion, in the Edit Boundary Condition window that appears, choose the option Fully fixed (U1=U2=U3=UR1=UR2=UR3=0).
   4. See Table 1 and Table 2^24,25,26 for all material property settings. Use the same material property, interaction relation, and boundary condition settings for other tissues and structures in pre- and postoperative models.
3. Validation of the individualized FE model
   1. Establish a loading point just posterior to the center of the upper endplate of the L3 vertebra before applying any loads. Couple all nodes on the L3 upper endplate with this loading point through constraint relationships.
   2. Apply different directions of pure bending moments of 3.5 N∙m at the model's loading point to simulate lumbar spine movements during flexion, extension, and lateral bending. Measure the ROM for each segment and compare it with the experimental data reported by Guan et al.²⁷.
   3. Apply a 150 N vertical load at the loading point and impose different directional loads of 2.5 N∙m, 5 N∙m, and 7.5 N∙m to simulate the motion of the lumbar spine in various directions. Measure the ROM for each segment and compare it with the experimental data reported by Panjabi et al.²⁸.
   4. Use the instantaneous axis of rotation method to measure and calculate the ROM for each lumbar segment.
      1. In the postprocessing module of the FE software, fix the view, capture before-and-after displacement images of the model in the same view, and import them into image processing software.
      2. Determine the instantaneous center and the spinal motion of each segment according to the methods described in the literature.
      3. For each segment, take measurements three times and use the average to minimize errors from the different measuring planes.
   5. Apply a 500 N vertical load and a 7.5 N∙m moment at the loading point to simulate flexion, extension, and lateral bending movements.
   6. In postprocessing, extract the maximum internal stress at the nucleus pulposus in the intervertebral discs of each segment and compare the data with the results reported by Dreischarf and Wike^14,29.

4. Loading of the FE model

1. Apply the same processes and load values to the preoperative and post-PLIF models to facilitate the analysis of the mechanical changes after PLIF surgery.
2. Apply a 400 N vertical downward load at the loading point above the L3 vertebra and a 7.5 N∙m moment load in each direction at the loading point to simulate human forward flexion, backward extension, lateral bending, and torsion motion.
   NOTE: Only motion in the unilateral direction during lateral bending and torsion needs to be simulated since the parametric lumbosacral model is symmetrical about the sagittal plane.
3. See Figure 1B for the final patient-specific post-PLIF FE model.

Representative Results

Simulation results of the patient-specific model compared to previous literature results
ROM of the intervertebral disc
According to the experimental loading conditions of Guan et al.²⁷, a pure bending moment load of 3.5 N∙m in different directions was applied at the loading point of the model to simulate the lumbar spine motion in flexion, extension, and lateral bending, and the ROM of each segment was measured and compared with the results of Guan's study. The comparison results are shown in Figure 2A-C. Compared with Guan's experimental data, the FE model established in this study had a smaller ROM for each segment in flexion and a larger ROM of L3-L4 in extension, both of which were basically within a reasonable range of experimental standard deviation. During lateral bending, the results of this study were all within the range of standard deviation.

For validation under a combined load (axial load and bending moment load), referring to the conditions of Panjiabi et al.²⁸in vitro experiment, a vertical load of 150 N was applied at the loading point, and moments of 2.5 N∙m, 5 N∙m, and 7.5 N∙m in different directions were applied to simulate lumbar spine motion in each direction. The results are shown in Figure 2D-I, where the direction of the y-axis indicates the direction of motion. Most of the results of this study are in good agreement with previous in vitro experimental data and the overall trend is similar. The ROM of the L4-L5 segment under a large load was slightly beyond the standard error range.

Stress in the nucleus pulposus
The results are shown in Figure 2J-L. We can conclude that the simulation results of the patient-specific model used here are aligned with those of other previously demonstrated validated FE models, as well as those of relevant in vitro experimental data^14,29.

Change in ROM of adjacent segments before and after PLIF surgery
The changes in the ROM of adjacent segments are shown in Figure 3. Significant increases in all directions of motion were found after fusion surgery. Intervertebral mobility of the L3-L4 and L5-S1 segments under forward flexion increased by 15.9% and 25.9%, respectively. Under posterior extension, the intervertebral mobility of the L3-L4 and L5-S1 segments increased by 5.9% and 15.6%, respectively. Under lateral bending, the mobility of the L3-L4 and L5-S1 segments increased by 10% and 17.5%, respectively, while the mobility of the L3-L4 and L5-S1 segments increased by 19% and 21.4%, respectively, during torsion, which is in line with the results of existing FE studies^30,31 and relevant in vitro experiments^32,33.

Changes in the overall ROM of the lumbosacral model before and after PLIF surgery
The overall lumbar mobility after L4-L5 fusion decreased by 33.5% during anterior flexion, 44.3% during posterior extension, 35.6% during lateral flexion, and 28.6% during torsion. It can be concluded that the overall ROM of the lumbosacral model decreased significantly in all directions of motion after fusion surgery, as shown in Table 3, indicating that although the mobility of the adjacent segments increased after fusion, the significant decrease in mobility of the fused segment resulted in a decrease in the overall ROM and an increase in the overall stiffness of the lumbosacral region.

Stresses in the facet joints of adjacent segments before and after PLIF
The average von Mises stress values of nine evenly distributed points of the facet joints were calculated on both sides; among them, the point with the highest value was selected as the index to compare the biomechanics of the facet joints before and after PLIF surgery.

As shown in Figure 4 and Figure 5, significant increases were found in the mean stresses in the facet joints of the adjacent segments after PLIF in all directions of motion. The mean stresses in the facet joints of the L3-L4 and L5-S1 segments increased by 42.2% and 45.3% during forward flexion, by 3.1% and 26.8% during extension, by 24.8% and 43% during lateral bending, and by 136.4% and 113% during torsion, respectively. These results are consistent with those in the literature^34,35. In this study, we also found that the increase in stress in the facet joints of the segments adjacent to the L3-L4 segment was slight during posterior extension (less than 5%), while the increase during torsion was extremely significant (an increase of over 100% in all instances).

Figure 5 shows the stress distribution in the facet joints of adjacent segments during motion in every direction before and after the PLIF procedure. The stresses in the facet joints of the L4-L5 segment decreased significantly after PLIF surgery, while the areas of highly (red area) and moderately (yellow and green areas) concentrated stress increased, indicating that not only did the maximum stress in the facet joints of the adjacent segments increase after fusion but also that the areas of concentrated stress increased, which might be one of the important factors influencing ASD.

Stress in the intervertebral discs of adjacent segments before and after PLIF
After PLIF surgery, the maximum stress in the annulus fibrosus of the L3-L4 and L5-S1 segments increased by 11.9% and 11.1% during forward flexion, 3.7% and 18.3% during extension, 47.6% and 59.5% during lateral bending, and 81.0% and 63.8% during torsion, respectively. The maximum stress in the nucleus pulposus of the L3-L4 and L5-S1 segments increased by 10.3% and 8.3% during forward flexion, 5% and 10.7% during extension, 32.3% and 21.6% during lateral bending, and 55.6% and 50% during torsion, respectively.

As shown in Figure 6, stresses in both the annulus fibrosus and the nucleus pulposus increased during motion in all directions after PLIF surgery, and the most significant increase occurred during torsion. Stresses in both the annulus fibrosus and the nucleus pulposus were greater at the L3-L4 segment than at the L5-S1 segment during motion in almost all directions preoperatively. After PLIF surgery, the most significant increase in stress occurred at the L3-L4 segment during torsion, with an 81% increase in the disc annulus fibrosus and 55.6% in the nucleus pulposus. The stress increase in the annulus fibrosus and nucleus pulposus during lateral bending was more significant than that during forward flexion and extension.

Figure 7 shows the stress distribution in the annulus fibrosus and nucleus pulposus of the adjacent segments before and after PLIF surgery. Similar to the pattern in the facet joints, stress concentration, and a more intensive stress distribution occurred in the annulus fibrosus and nucleus pulposus during torsion. In addition, the maximum internal stress during anterior flexion, extension, and lateral bending movements was found on the side being compressed, whereas the location of the maximum stress during torsion was not confined to one side.

Figure 1: Parameters and the procedure for developing patient-specific FE models. (A) Parameters for generating patient-specific lumbar spine geometry. (B) Procedure for developing patient-specific preoperative and postoperative FE models of the lumbosacral spine (L3-S1). Please click here to view a larger version of this figure.

Figure 2: Simulation results of the patient-specific model compared to previous literature results. (A-C) ROM of the patient-specific model compared to that of Guan's study. (D-F) Comparison of ROM between FE analysis using the patient-specific model and Panjiabi's data under a similar combined load at the L3-L4 segment. (G-I) Comparison of ROM between FE analysis using the patient-specific model and Panjiabi's data under a similar combined load at the L4-L5 segment. (J-L) Results of stress simulation in the nucleus pulposus between the present study and other studies. Please click here to view a larger version of this figure.

Figure 3: ROM of the L3-L4 segment and the L5-S1 segment before and after PLIF during motion in different directions. (A) L3-L4 segment and (B) L5-S1 segment. Please click here to view a larger version of this figure.

Figure 4: Mean stress during flexion, extension, lateral bending, and torsion in facet joints before and after PLIF. (A) Flexion. (B) Extension. (C) Lateral bending. (D) Torsion. Please click here to view a larger version of this figure.

Figure 5: Stress distribution in facet joints during motion in different directions before and after PLIF. Please click here to view a larger version of this figure.

Figure 6: Internal stress in the annulus fibrosus and nucleus pulposus of adjacent segments before and after PLIF. (A) Annulus fibrosus. (B) Nucleus pulposus. Please click here to view a larger version of this figure.

Figure 7: Stress distribution in the annulus fibrosus and nucleus pulposus of adjacent segments before and after PLIF surgery. Please click here to view a larger version of this figure.

Instrument	Material properties	Mesh type	Young's modulus (MPa)	Poisson's ratio
Pedicle screw	Titanium	C3D4	110000	0.3
Fixation rods
Intervertebral fusion device	Polyetheretherketone	C3D4	3700	0.3

Table 1: Material properties of internal fixation and fusion instrumentation

Structure	Unit type	Modulus of elasticity (MPa)		Poisson's ratio	Density (kg/mm3)
Cortical bone	S4	12000		0.3	1.7 × 10-6
Cortical bone	C3D4	100		0.2	1.1 × 10-6
Posterior element	C3D4	3500		0.25	1.4 × 10-6
End plate	S4	23.8		0.4	1.2 × 10-6
Annulus fibrosus	C3D8H	C10 = 0.18, C01 = 0.045		-	1.05 × 10-6
Nucleus pulposus	C3D8H	C10 = 0.12, C01 = 0.03		-	1.02 × 10-6

Table 2: Material property parameters of the lumbar FE model

Moving directions	Maximum displacement before PLIF/mm	Maximum displacement after PLIF/mm	Percentage change
Forward flexion	49.8	33.1	-33.50%
Extension	19.2	10.7	-44.30%
Lateral bending	29.8	19.2	-35.60%
Torsion	18.5	13.2	-28.60%

Table 3: Maximum displacement during motion in different directions before and after PLIF

Discussion

In this study, a geometrically parametric patient-specific FE model was established to analyze the biomechanical characteristics of the lumbar spine after PLIF surgery. The results showed that the stress in the facet joints and disc of the fused segment decreased significantly after PLIF surgery, indicating that PLIF could effectively strengthen the stability of the decompressed segment and delay further aggravation of the lesion. The overall mobility of the lumbar spine decreased after PLIF surgery, while the ROM, facet joint stress, areas of concentrated stress in the facet joints, and stress in adjacent intervertebral discs all increased to varying degrees. This suggests that fusion not only restricted lumbar motion but also allowed the adjacent segments to move over a larger range and share more of the load to compensate for the fused segment. A vicious cycle, therefore, forms that further increases the stress in the adjacent segments, resulting in increasing stress concentration and increasing risk of degenerative disease in the adjacent segments. Our results are consistent with those of FE studies by other researchers^36,37, and cadaver specimen experiments by Weinhoffer²⁶ and Cunningham³⁸. Although the magnitude of stress increase in the intervertebral disc varies widely from study to study, we believe this variation is related to differences in the samples, loading, or measurement methods.

We also found that stresses in the facet joints of the proximal adjacent segment (L3-L4) were higher than those in the facet joints of the distal segment (L5-S1) during motion in almost all directions (forward flexion, backward extension, and torsion). The stress increase in both the proximal and distal adjacent segments during torsion was more significant (>50%) than that during motion in other directions. In particular, the stress in the facet joints and annulus fibrosus of the proximal adjacent segment (L3-L4) increased to a larger extent than that in the distal adjacent segment (L5-S1). It has been suggested that the proximal adjacent segment is more prone to degeneration than the distal adjacent segment after spinal fusion due to the significant increase in stress³⁹.

Based on the findings of this study, it is critically important for surgeons to implement strategies during spinal fusion surgeries to reduce the risk of ASD, particularly in the proximal segment. This includes preoperative planning, such as assessing risk factors for ASD, e.g., pre-existing degeneration in adjacent segments or abnormal sagittal balance. Additionally, adopting a more conservative approach during the surgical process, such as limiting the number of fused segments, is advisable. Moreover, considering the significant increase in stress during torsion, proactive measures such as fusing the facet joints of the superior adjacent segment or encouraging patients to compensate for lumbar motion with lower limb activities during postoperative recovery to minimize the amplitude of lumbar torsion can be beneficial.

This study proposes a method for individualized parameterization in lumbosacral FE modeling. The validation test performed in this study demonstrated that the established L3-S1 lumbosacral FE model in this study aligns with the data in previous literature^27,28. A small amount of simulation data slightly exceeded the standard error range, such as the ROM of the L4-L5 segment under a large load, which may have been caused by individual differences or model parameter differences. However, the deviation was not large, and the trend was in line with the results of the literature. Moreover, we believe that this difference has little impact on the simulation results in other directions and the final clinical explanation, given that the human body rarely bears large loads under horizontal torsion. Regarding stress in the nucleus pulposus under pressure, the simulation results of the model used in this study were also aligned with the literature data^14,29. Generally, the biomechanical characteristics of this study model are consistent with the existing theory, which indicates that the model can be used to simulate the motion of the lumbosacral region under conventional loads as well as to analyze the individual biochemical characteristics after lumbar spinal surgery.

The crux of this modeling lies in the use of simple geometric shapes to closely mimic the form of the vertebrae and surrounding tissues, particularly the shape of the load-bearing bones. This is achieved in the modeling software by creating a parameter-controlled model composed of simple geometries through straightforward operational methods. The main challenge in this process is positioning the modeled spinal segments in three-dimensional space to form a coherent lumbosacral spinal model using lateral X-ray images as a sagittal reference. This requires manual adjustment on the sagittal, coronal, and transverse planes to avoid spatial misalignment or interference between vertebrae. Currently, this step can only be completed manually, which is time-consuming and cannot be easily resolved without the introduction of feature point recognition technology.

The modeling method proposed in this study also has certain limitations. First, this parameterized model primarily considers internal fixation or fusion in the lumbar spine, with significant simplifications made to the morphology of the sacrum. This is because, in studies focusing on the lumbar spine's load-bearing characteristics, the sacrum is typically assumed to be immobile. Hence, the model is not applicable in situations with sacral involvement in internal fixation. Second, each vertebral body in the model is separately defined with cortical and cancellous bone, assuming uniform thickness and homogeneity for the cortical bone and homogeneity for the cancellous bone. This does not account for the potential variations in bone density caused by osteoporosis. Third, while this study employs a simplified static modeling approach, predominantly using basic geometric shapes for representation of the lumbosacral region, there are more complex models^20,40, such as viscoelastic, poroelastic, and poro-hyperelastic models commonly found in existing lumbar spine studies. The choice of a simplified static model was driven by our focus on individualized geometric shapes within lumbosacral FE modeling. This approach, while less complex than that of viscoelastic or poroelastic models, offers significant advantages in terms of computational efficiency and ease of individual customization. It allows for rapid adaptation to individual anatomical variations, which is critical for personalized clinical applications. Compared to more complex viscoelastic or poroelastic models, the simplified static model used here may not capture certain dynamic or fluidic aspects of spinal biomechanics with the same fidelity. However, for the purposes of this study - focusing on individualized geometric shapes and conventional load simulations - the model provides a balanced trade-off between detail and practical applicability.

In this study, the geometrically parametric FE modeling method showed many advantages, such as less computation time, better model convergence, and more convenience in achieving patient-specific modeling than the CT modeling method used at present. Future studies could focus on merging the modeling method used here with image recognition based on artificial intelligence algorithms to amplify fully automated modeling and efficient biomechanical analysis. Our modeling method could also be applied in biomechanical analysis studies with a large sample size, thus improving the reliability of the conclusions while elucidating the mechanisms of spinal diseases to develop innovative treatment methods and prevention strategies.

Materials

Name	Company	Catalog Number	Comments
Abaqus	Dassault	https://www.3ds.com/products/simulia/abaqus	Finite element analysis
AutoCAD	Autodesk	https://www.autodesk.com/products/autocad/	An Engineering Computer Aided Design software used to measure the ROM of different vertebral segment
CT scan dataset	China Japan Friendship Hospital		Dataset of an adult healthy male with no history of trauma, deformity or tumor of the spine (height 180 cm, weight 68 kg).The raw data were stored in Dicom 3.0 format with a pixel size of 0.33 mm and a layer spacing of 1 mm.
Hypermesh 2019	Altair	https://altair.com/hypermesh/	Mesh generation
Mimics Research 21.0	Materialise	https://www.materialise.com/en/healthcare/mimics-innovation-suite/mimics	Model construction