Institutional review board (IRB) approval for the imaging procedures used in this study was obtained from the IRB of Medical College of Wisconsin, Milwaukee, USA (IRB approval no., PRO00043133).

This study utilized preoperative MRI from a volunteer diagnosed with mild DCM, with a modified Japanese Orthopedic Association score of 16 (Fig. 1). The MRI data was acquired using a 3T GE Premiere machine (GE Healthcare, Chicago, IL, USA) capturing cervical spine T2-weighted fast spin echo images with the following parameters: TR/TE (repetition time/echo time) of 2,500/122 ms and a 3D isotropic resolution of 0.8×0.8×0.8 mm. The spinal geometry was extracted from the MRI using the RadiAnt DICOM viewer software (Medixant, Poznan, Poland). Metrics such as the dimensions of vertebral height, width, and depth, and the length of the spinous process and lateral mass width, were measured linearly. Additional metrics included the vertebral body-to-spinous process length, disk height, Cobb angle, and the thickness of the anterior and posterior longitudinal ligaments. The degeneration of each cervical disk was categorized based on the T2-weighted MRI images [14]. Patient-specific geometries were integrated into a previously validated generic Finite element (FE) model, creating a patient-specific cervical spine FE model [15,16].

The FE model of the C2–T1 cervical spine was developed using published material properties (Table 1), as previously described [6,17,18]. This model represented a 60-year-old man with cervical myelopathy, characterized by spinal cord compression due to a disk-osteophyte complex at C4–5, C5–6, and C6–7. Bilateral foraminal stenosis was present at each level, with maximum cord compression at C5–6, accompanied by T2-hyperintensity within the spinal cord. Since bone mineral density values for the cervical vertebrae were unavailable, literature-derived bone density values were employed [19,20]. For the vertebral structure, we employed hexahedral elements constructed from isotropic linear elastic materials to represent cancellous bone. Conversely the cortical shell, a 0.5 mm thick layer, was modeled using quadrilateral elements and isotropic linear elastic materials, mirroring the vertebral body’s geometry. The 0.2 mm thick endplates utilized the same element type and material properties as the cortical shell. The intervertebral disk was segmented into three components: nucleus, annulus ground, and annulus fibers. Connections between nodes were established for both the disk and endplates, as well as for the cancellous bone and endplates. Annulus fibrosus was modeled using a nonlinear, orthotropic material model, with fibers arranged in distinct patterns depending on their anterior or posterior location. Our model also featured five primary ligament groups: the anterior longitudinal ligament, posterior longitudinal ligament, interspinous ligament, ligamentum flavum (LF), and capsular ligament (CL). These ligaments were modeled using quadrilateral membrane elements with nonlinear stress-strain relationships to enhance realism (Table 1). The material models for the annulus ground and nucleus pulposus were designated as hill foam and viscoelastic fluid, respectively. These components were meshed using hexahedral and quadrilateral elements, with a total of 44,799 hexahedral and 35,679 quadrilateral elements employed in the overall model. This morphable generic FE model was validated against data from human cadaver cervical columns. This validation confirmed that the predicted range of motion (ROM) at all segments fell within one standard deviation of the mean observed in cadaveric samples, ensuring the model’s accuracy [21,22].

Spinal cord geometry was assessed using T2 MRI. Geometric measurements were taken at anatomical landmarks such as the superior endplate, inferior endplate, disk segment, and mid-vertebral sections. Using RadiAnt software (Medixant, Poznań, Poland), the two-dimensional geometry of the spinal cord and surrounding cerebrospinal fluid (CSF) were outlined. These geometrical outlines were then transferred into CATIA V6 (Dassault Systemes, Vélizy-Villacoublay, France) for solid model creation, employing spline tools and interpolation techniques for accurate representation.

For meshing, a hexahedral grid was generated using the ANSA software (BETA CAE Systems, Farmington Hills, MI, USA). Our spinal cord model integrated various anatomical components, including gray and white mater, as well as surrounding tissues. Specifically, the pia mater covered the outer surface of the white mater, while the dura mater lined the outer layer of the CSF. Filling the gap between the pia mater and the dura mater was CSF. Denticulate ligaments were modeled according to existing anatomical literature and were attached to the lateral surface of the pia mater. Node-to-node connections ensured seamless integration across all anatomical features. The patient-specific spinal cord FE model was then integrated into the patient-specific FE model of the osteoligamentous cervical spine. Within the spinal canal, the spinal cord was left unconstrained. Frictionless surface-to-surface interactions were set between the spinal column and spinal cord to simulate physiological conditions. Using this patient-specific FE model, we simulated two surgical interventions ACDF (C4–7) and PCLF (C5–6 laminectomy and C4–C7 fusion).

For the ACDF procedure, a central bone graft (12.5 mm×15 mm), was inserted between the C4–C5, C5–C6, and C6–C7 vertebral bodies. The graft was modeled using trabecular bone material properties. To ensure precise representation, 3D models of the anterior cervical plate and variable angle screws were created using the CATIA V6 software (Dassault Systemes). The cervical plate measured 58.2 mm in height, 17.5 mm in width, and 2 mm in thickness, spanning from C4 to C7 along the anterior cortex of the vertebrae. Standardized screws (18 mm in length and 3 mm in diameter) were simulated in the model. The top and bottom screws were inserted at a 35° angle for optimal fixation, while the remaining screws were inserted parallel to the vertebral endplates. Tied constraints were applied to simulate the integration of the bone graft with adjacent vertebral bodies. Identical constraints modeled the interaction between screws and vertebral bodies [23]. Automatic surface-to-surface constraints definition was used to model the interaction between screws and titanium plate, as well as between the plate and vertebral bodies. To enhance clinical relevance, the model also simulated spinal cord decompression after ACDF. This included creating a ventral CSF column of approximately 1 mm at the decompressed segments, replicating the postoperative anatomical changes.

PCLF was modeled by removing the spinous process, lamina, and the LF, along with the interspinous ligaments at C5 and C6 [12]. The interspinous ligaments and ligamentous flavum were also removed at the adjacent segments (C4–C5 and C6–C7). Subsequently, fusion from C4 to C7 was modeled using the CATIA V6 software. For the PCLF surgery, custom-designed titanium implants were created, featuring a lateral mass screw with a 3.5 mm diameter and 14 mm length, and a connecting rod with a 3.5 mm diameter and curvature ranging from 0° to 20°. Set screws and screw heads were designed to secure the rod to the vertebrae. Screws were inserted parallel to the facet joints from C4–C7 using mono-cortical fixation. The rod was manipulated to ensure optimal engagement with the lateral mass screws. Constraints were applied to simulate a rigid union between the screws and the rod while an automatic surface-to-surface contact definition was used to model the interaction between the screws and the plate. To simulate the biomechanical changes resulting from PCLF surgeries, the spinal cord was repositioned dorsally by an average of 2 mm at the decompressed segments [24]. This mimicked the expected dorsal spinal cord following posterior decompression performed for DCM.

Biomechanical responses of the spine and spinal cord were analyzed under various loading conditions, including flexion, extension and lateral bending, and axial rotation using LS-PrePost ver. 4.3 and LS-Dyna ver. 9.0.1 (ANSYS Inc., Livermore, CA, USA). The models were constrained at the inferior surface of the T1 vertebra in all six degrees of freedom. A pure moment of 2 Nm was applied in both flexion and extension modes and for lateral bending and axial rotation [21]. A follower load of 75 N was also applied during flexion and extension to simulate head mass and muscle force [25]. This follower load was executed through strategically placed cables, with node locations guided by each segmental unit’s center of rotation. No follower load was applied during lateral bending and axial rotation, as recommended in prior studies [22]. Outcomes included segmental ROM, spinal cord stress (von Mises Stress), and strain (Maximum Principal Strain). These metrics were evaluated at both the fused and adjacent spinal segments. For each surgical simulation, additional metrics such as CL strain and intradiskal pressure were evaluated and compared with the preoperative patient-specific FE model. All recorded average stress/strain values were filtered to exclude data below the 5th percentile and above the 95th percentile [26].

The intact C2–T1 FE model was validated by comparisons with angular motions derived from human cadaver experiments, focusing on segments C2–C3 to C7–T1. Validation results demonstrated that the FE model’s response to applied loads fell within ±1 standard deviation of the ROM values from experimental data, demonstrating the model’s accuracy [21,22].

The biomechanical responses of the ACDF model were compared to the cadaveric study by Tsitsopoulos et al. [27], while the PCLF model was validated against the findings of Kretzer et al. [28]. Biomechanical responses derived from both FEMs were within 1 standard deviation of the published cadaver data (Fig. 2).

Segmental ROM was reduced at the fused segments in both ACDF and PCLF models, with no substantial difference observed under physiological loading conditions (Fig. 3). Preoperative ROM averages across all segments ranged from 2.2° to 3.1° for flexion, extension, lateral bending, and axial rotation. In both the ACDF and PCLF models, there was a small reduction in the average ROM, with values ranging from 1.6° to 2.7° and 1.7° to 2.8°, respectively. At the superior adjacent segment, the ACDF model exhibited a higher ROM during flexion compared to PCLF (32.8% versus 7.7%). Conversely, in extension, PCLF showed a higher ROM than ACDF (22.9% versus 1.8%). At the inferior adjacent segment, both models displayed similar increases in ROM during flexion (18.6% versus 17.7%). However, during extension, PCLF demonstrated a higher ROM than ACDF (18.2% versus 13.3%).

Fig. 4 illustrates the changes in spinal cord stress and strain at the fused segments (C4–C7). Notable differences between ACDF and PCLF models were observed. At the C4–C7 segments, ACDF exhibited a greater reduction in spinal cord stress during flexion compared to PCLF (−40% versus −13%). In contrast, during extension, PCLF showed a greater stress reduction than ACDF (−61% versus −28%). Regarding spinal cord strain at C4–C7 segments, both models demonstrated similar reductions during flexion (−18% versus −20%), while PCLF achieved a more substantial decrease during extension compared to ACDF (−57% versus −23%).

At the superior adjacent segment, the ACDF model showed greater spinal cord stress than PCLF during flexion (67.1% versus 10.5%). In contrast, during extension, PCLF exhibited a higher stress than ACDF (39% versus −6.8%). For spinal cord strain at this segment, ACDF showed higher strain during flexion (24.1% versus 10.3%) compared to PCLF, while PCLF exhibited more strain during extension (45.4% versus 9%). At the inferior adjacent segment, ACDF displayed higher spinal cord stress than PCLF during flexion (43.7% versus −3.5%). During extension, both models showed decreased stress (−16.4% versus −26.2%). For spinal cord strain, both models exhibited slight decreases during flexion (−3% versus −1.8%), while PCLF demonstrated a 19.2% increase during extension, with no change observed in ACDF.

At the superior adjacent segment, the ACDF model showed higher intradiscal pressure during flexion compared to PCLF (35.1% versus 28.1%), while PCLF exhibited higher intradiscal pressure during extension compared to ACDF (33.5% versus 16.3%). At the inferior adjacent segment, ACDF displayed a higher intradiskal pressure in flexion than PCLF (23.2% versus 16.4%), while in extension, both models exhibited similar pressures (26.4% versus 25.7%). Fig. 5 displays the distribution of intradiskal pressure across the annulus and nucleus.

At the fused segments, both models demonstrated a reduction in average CL strain, with PCLF showing greater reductions in axial rotation strain compared to ACDF (−84.2% versus −69.4%). At the superior adjacent segment, ACDF exhibited higher strain than PCLF (37.9% versus 6.9%) during flexion, whereas in extension, PCLF had a higher strain compared to ACDF (24.8% versus −0.2%) (Fig. 3).