A three-dimensional (3D) nonlinear FE model of the lumbar spine was developed which consisted of three lumbar vertebrae, two intervertebral discs, and the associated spinal ligaments. Geometrical details of the human lumbar spine (L2-4) were obtained from high-resolution computed tomography (CT) images of a 46-year-old male subject who had no spinal deformities. Digital CT data were imported into a software program (Mimics; Materialise Inc., Leuven, Belgium) and used to generate a 3D geometrical surface of the lumbar spine. Exported IGES files from the Mimics software were input into Unigraphics NX 3.0 (Siemens PLM Software, Torrance, CA, USA) to form solid models for each vertebral segment. The solid model was then imported into Hypermesh 8.0 (Altair Engineering, Inc., Troy, MI, USA) to generate FE meshes. The FE method was analyzed with commercially available software (ABAQUS 6.6-1; Hibbitt, Karlsson and Sorenson, Inc., Providence, RI, USA).16,19

3D homogenous and transversely isotropic solid elements were used to model the cortical and cancellous cores and the posterior bony parts of the vertebrae. The anterior longitudinal ligament, posterior longitudinal ligament, intertransverse ligament, ligament flavum, capsular ligament, interspinous ligament, and supraspinous ligament were modeled using tension-only truss elements.16,19

Material properties were selected based on various literature sources (Table 1).20,21,22,23,24 The cortical and cancellous regions of the vertebrae were modeled independently. Differentiating between cortical and trabecular bone in the posterior region was difficult; therefore, the posterior elements were all assigned a single set of material properties.

The annulus fibrosus (AF) was modeled as a composite of a solid matrix with embedded fibers (via the REBAR parameter) in concentric rings surrounding a nucleus pulposus (NP), which was considered to be an incompressible inviscid fluid. Element members with hybrid formulation (C3D8H) combined with a low elastic modulus and large Poisson ratio definitions were applied to simulate the NP. Eight-node brick elements were employed to model the matrix of the ground substance. Each of the four concentric rings of the ground substance contained two evenly spaced layers of annulus fibers oriented at ±30° to the horizontal plane. The reinforcement structure annulus fibers were represented by truss elements with modified tension-only elasticity. In the radial direction, four double cross-linked fiber layers were defined, and those fibers were bounded by the annulus ground substance and both endplates. In addition, these fibers proportionally decreased the elastic strength from the outermost (550 MPa) to the innermost (358 MPa) layer.16,19,25,26

The articulating facet joint surfaces were modeled using surface-to-surface contact elements in combination with the penalty algorithm with a normal contact stiffness of 200 N/mm and a friction coefficient of zero. The thickness of the cartilage layer of the facet joint was assumed to be 0.2 mm. The initial gap between the cartilage layers was assumed to be 0.5 mm. The cartilage was assumed to be isotropic: a linear elastic with a Young's modulus of 35 MPa and a Poisson's ratio of 0.4.27 Spinal ligaments were represented with nonlinear material properties. Naturally-occurring changes in ligament stiffness (initially low stiffness at low strains, followed by increasing stiffness at higher strains) were simulated through a "hypoelastic" material designation (Table 1). 3D truss elements were used to simulate ligaments, which were active only upon tension.16,19

Two kinds of decompression states were simulated. One was hemi-laminectomy at the L3-4 segment, for which the right lamina was totally removed at the isthmic portion. The other was total (bilateral) laminectomy. In order to simulate the total laminectomy, both sides of the lamina were removed at the isthmic portion. In both models, a supraspinous ligament and interspinous ligament between the L3 and L4 spinous processes were removed along with a part of the L3 and L4 spinous processes, and the ligamentum flavum of L3-4 was also removed depending on the extent of the laminectomy. All models had continuity of the proximal posterior ligament complex between the L2 spinous process and the remaining L3 spinous process (Fig. 1).

A posterior pedicle screw fixation was simulated in the L3-4 intertransverse process fusion. All screws had sharp threads to prevent relative motion at the bone-screw interface and were inserted into the anterior third of the vertebral body through the pedicles. The pedicle screws did not violate the medial wall of the pedicle or the endplates. With the exception of the screw tip, the remaining surface of the screw was fixed to the bone without allowing relative motion. A "tie" contact condition was used to enable the screw threads and vertebrae to be permanently bonded together by full constraint. The diameters of all pedicle screws were assumed to be 5.0 mm, with a mean outer diameter of 6.5 mm (including thread height). The lengths of the screws in the L3 and L4 were 40 mm (Fig. 2). For the simulation of unilateral pedicle screw fixation, the pedicle screw-rod system was applied only at the right side, while the pedicle screw-rod system was applied to both sides for the bilateral pedicle screw fixation simulation.16

Four scenarios were simulated in total according to the extent of decompression or pedicle screw fixation. The BF-HL model represented hemi-laminectomy with bilateral pedicle screw fixation in the L3-4 segment, and the BF-TL model represented total laminectomy with bilateral pedicle screw fixation. Furthermore, in the UF-HL model, hemi-laminectomy was simulated with unilateral pedicle screw fixation and total laminectomy was simulated with unilateral pedicle screw fixation in the UF-TL model (Fig. 1).

This FE investigation included two types of loading conditions corresponding to 1) loads used in the experimental part of the studies28,29,30 for model validation and 2) model predictions for clinically relevant loading scenarios. For model validation, the experimental and simulated loading protocols were identical. For four pure moments, the results of range of motion were compared to those of previous in vitro studies by Renner, et al.28 and Wilke, et al.30 The nodes of the inferior surfaces of the most inferior vertebral body were completely fixed in all directions. To validate the model for pure moments, we applied the same loading conditions as Renner, et al.28 and Wilke, et al.,30 respectively. Nodes on top of the L2 vertebra were defined as coupling nodes. A reference node was created and connected to all coupling nodes. A coupling element was created to distribute moments on the reference node. For comparison with the results of study by Renner, et al.,28 8 Nm in flexion and 6 Nm in extension were imposed on the L2 vertebral body without preload and with an 800 N follower preload. ±4 Nm in torsion and ±6 Nm in lateral bending moments were imposed on the L2 vertebral body without preload. For comparison with the results of study by Wilke, et al.,30 7.5 Nm were imposed on the L2 vertebral body in flexion, extension, lateral bending, and torsion moment, respectively. Furthermore, the model was also validated for using intradiscal pressure. The intradiscal pressure of the current model was compared with the results in the experimental study by Schilling, et al.29 The experimental and simulated loading protocols were identical.

For model prediction for clinically relevant scenarios, the loading condition followed the hybrid testing protocol implemented during the flexibility testing of the FE models as previously described for the study of adjacent level biomechanics.31 The follower load technique was used to simulate the vector sum of the trunk muscle co-activation using a single internal force vector acting tangent to the curvature of the spine passing through each segmental center of rotation.32 The 400 N compressive follower load was simulated at each motion segment in the model by a pair of 2-node thermo-isotropic truss elements. The trusses were attached bilaterally to the cortical shell of the vertebrae at each motion segment. Each truss spanned the disc space passing through the instantaneous center of rotation at each motion segment.28 The hybrid testing protocol involved the application of the pure moment to the intact and fusion models until the L2-4 rotation (displacement) equaled the intact load control case values, which was achieved by imposing 7.5 Nm flexion, extension, lateral bending and torsion moments on the L2 vertebral body in the intact model.

For model validation, the experimental and simulated loading protocols were identical. For four pure moments, the results of range of motion were compared to those of previous in vitro studies by Renner, et al.28 and Wilke, et al.30 (Fig. 2A and B). For the intradiscal pressure, the results were compared to those in the study by Schilling, et al.29 (Fig. 2C). We performed a mesh convergence test with spinal segment (L2-4). As critical result parameters, mesh densities about range of motion were selected to test for difference among the intact models. The finest mesh density was chosen since the change was within 1.8% (<0.2°). Finally, the intact spine model was converged using 52955 nodes and 47213 elements.

The range of motion (ROM) at each corresponding level was compared among the 4 fusion models. The change of ROM at the adjacent segment was described as the percent change from the intact model in each moment under a hybrid protocol, and the change of ROM at the fusion segment was described as the percent of remaining ROM (100-percent change) from the intact model in each moment under a hybrid protocol (Figs. 3 and 4).

Under flexion moment, the ROM at the adjacent segment increased by 15.4%, 15.4%, 7.69%, 7.69% in the BF-HL, BF-TL, UF-HL, and UF-TL models, respectively, compared to the intact model (Fig. 3). Under extension moment, the ROM at the adjacent segment increased by 16.7%, 16.7%, 8.3%, 8.3% in the BF-HL, BF-TL, UF-HL, and UF-TL models, respectively, compared to the intact model. Therefore, under flexion and extension moments, the unilateral fixation leads to about a 50% reduction in the ROM increase at the adjacent segment after fusion when compared to bilateral fixation. Under torsional moments, the ROM at the L2-3 segment increased by 14.9%, 14.9%, 9.9%, and 4.9% in the BF-HL, BF-TL, UF-HL, and UF-TL models, respectively, compared to the intact model. Under lateral bending moments, the ROM at the adjacent segment increased by 6.25%, 6.25%, 4.17%, and 2.08% in the BF-HL, BF-TL, UF-HL, and UF-TL models, respectively, compared to the intact model (Fig. 3). Therefore, the extent of decompression influenced the percent increment only in the unilateral fixation models and only under torsion and lateral bending.

Relatively larger motions were noted at the fusion segment (L3-4) in the unilateral fixation (UF-HL and UF-TL) models (Fig. 4). Under flexion moments, the remaining range of motion was 8.9%, 8.9%, 32.1%, and 32.1% in the BF-HL, BF-TL, UF-HL, and UF-TL models, respectively, compared to the intact model. Under extension moments, 16.0% motion remained in the BF-HL and BF-TL models, as compared to the intact model. However, 47.8% motion remained in the UF-HL and UF-TL models. Under torsion moments, bilateral fixation models (BF-HL and BF-TL) showed 50.5% remaining motion at the fusion segment, while 61.1% and 77.7% range of motion remained at the fusion segment in the UF-HL and UF-TL models, respectively. Under lateral bending, the remaining range of motion was 26.8%, 26.8%, 31.7%, and 34.1% in the BF-HL, BF-TL, UF-HL, and UF-TL models, respectively, compared to the intact model (Fig. 4).

The four models showed increased maximal von Mises stress on the AF at superior adjacent L2-3 segments under the four moments, similar to the increases in ROM. The bilateral pedicle screw models, BF-HL and BF-TL, experienced the largest increase in AF stress at the adjacent L2-3 segments under all four moments. Compared to the intact model, both the BF-HL and BF-TL models yielded increases in AF stress of 34.1%, 125.0%, 84.2%, and 87.5% at the L2-3 segment under moments of flexion, extension, torsion, and lateral bending, respectively (Fig. 5). Following the bilateral pedicle screw fixation models, both the UF-HL and UF-TL models demonstrated 23.1% and 112.5% increases in AF stress under flexion and extension moments, respectively. Furthermore, under torsion and lateral bending moments, the unilateral fixation models demonstrated different increments in AF stress related to the extent of decompression. The UF-HL and the UF-TL model demonstrated 71.2% and 57.9% increases in AF stress, respectively, under torsion moments. Under lateral moments, the UF-HL and the UF-TL model showed 50.0% and 37.5% increases, respectively. Generally, unilateral fixation models led to a prominent reduction of increased AF stress at the adjacent L2-3 segment under the four moments (Fig. 5).