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

Summary

Here, we used a patient-specific finite element model to analyze the mechanical changes in adjacent segments after spinal fusion surgery. The results showed that fusion surgery reduced the overall motion of the lumbar spine but increased the load on and stress in adjacent segments, especially the proximal segment.

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.

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 G...

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...

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