The establishment of intact lumbar spine model (INT model)

To create a three-dimensional FE model, computed tomography scan DICOM files of the L1 to L5 lumbar spine of a middle-aged male were obtained at 1-mm intervals. The commercially available visualization software Amira 3.1.1 (Mercury Computer Systems, Inc., Berlin, Germany) was used to describe cross-section contours of each spinal component in accordance with gray scale value (Figure 3.1). Then, the three-dimensional surface geometries were constructed through sequential processed cross-section contours as shown in Figure 3.2 A.

Each spinal component was exported as a Drawing Exchange Format (DXF) file and converted to the Initial Graphics Exchange Specification (IGES) file as shown in Figure 3.2 B. The FE analysis software ANSYS 14.0 (ANSYS Inc., Canonsburg, PA) was used to reconstruct the FE model by converting the IGES file to ANSYS Parametric Design Language (APDL) code in Figure 3.2 C. The INT model was an osseo-ligamentous lumbar spine, which included the vertebrae, intervertebral discs, endplates, posterior bony elements, and all seven ligaments (Figure 3.3).

27

Eight-node solid element (SOLID185) were used for modeling the cortical bone, cancellous bone, posterior bony element, cartilage endplate, and annulus ground substance. The cortical bone and cancellous bone were assumed to be homogeneous and transversely isotropic [90]. The posterior bony element and cartilage endplate were assumed to be homogeneous and isotropic [90]. The intervertebral disc consisted of annulus ground substance, nucleus pulposus and collagen fibers embedded in the ground substance. The nonlinear annulus ground substance was simulated by using a hyper-elastic Mooney-Rivlin formulation [91, 92]. The collagen fibers simply connected between nodes on adjacent endplates to create an irregular criss-cross configuration. These irregular angles of collagen fibers were oriented within the range of the Marchand’s study[93]. In the radial direction, twelve double cross-linked fiber layers were defined to decrease elastic strength proportionally from the outermost layer to the innermost.

Therefore, the collagen fibers in different annulus layers were weighted (elastic modulus at the outermost layers 1-3: 1.0, layers 4-6: 0.9, layers 7-9: 0.75, and at the innermost layers 10-12:

0.65; cross sectional areas at the outermost layers 1-3: 1.0, layers 4-6: 0.78, layers 7-9: 0.62, and at the innermost layers 10-12: 0.47) based on previous studies [94, 95]. The nucleus pulposus was modeled as an incompressible fluid with a bulk modulus of 1666.7 MPa by eight-node fluid elements (FLUID80) [90]. 43 % of the cross-sectional area in the disc was defined as the nucleus, which was within the range of the study by Panagiotacopulos (30-50 %) [96]

Therefore, approximately 47 % to 49 % disc volume was assigned to nucleus pulposus. All seven ligaments and collagen fibers were simulated by using two-node bilinear link elements (LINK10) with uniaxial tension resistance only, which were arranged in an anatomically correct direction [97]. The cross-sectional area of each ligament was obtained from previous studies [94, 98, 100], and material properties of the spine are listed in Table 3.1. The facet joint was treated as having sliding contact behavior using three-dimensional eight-node surface-to-surface contact elements (CONTA174), which may slide between three-dimensional target elements (TARGE170). The coefficient of friction was set at 0.1 [101]. The initial gap

28

between a pair of facet surfaces was kept within 0.5 mm [90]. The stiffness of the spinal structure changes depending on the contact status, so the standard contact option in ANSYS was adopted to account for the changing-states nonlinear problem in this study. In addition, the element’s shape will change after applying bending moments, thus changing the individual element stiffness. Therefore, the large displacement analysis option in ANSYS was chosen to solve this geometric nonlinear problem. The INT model consisted of 112,174 elements and 94,162 nodes [102, 103].

Figure 3. 1: Each spinal component was selected from computed tomography scan DICOM file to create material-related contours.

29

Figure 3. 2: Modeling process of the L3 vertebra: (A) surface geometries of vertebra were reconstructed through sequential processed computed tomography scan DICOM file; (B)

surface geometry was exported to the DXF file; (C) FE model of the L3 vertebra.

Figure 3. 3: The lumbar finite element model used in this study (Intact model from L1 to L5 levels)

30

Table 3. 1: Material properties used in the FE model

Material Element type Young’s modulus

(MPa)

Poisson’s ratio Area (mm²)

References

Vertebral

Cortical 8node-Solid 185 Ex=11300

Ey=11300

Cancellous 8node-Solid 185 Ex=140

Ey=140

Annulus fibers 2node-Link 10 [94, 95]

Outmost (1-3 layers) 550 - 0.76

*ALL, anterior longitudinal ligament; PLL, posterior longitudinal ligament; TL, transverse ligament; LF, ligamentum flavum; ISL, interspinous ligament; SSL, supraspinous ligament; CL, capsular ligament.

3.1.1.1. Convergence test of INT model

In order to get reliable data, convergence test were conducted. Three mesh densities (coarse model: 4,750 elements / 4,960 nodes; normal model: 27,244 elements / 30,630 nodes;

31

finest model: 112,174 elements / 94,162 nodes) were selected to test ROM in the intact model (Figure 3.4). The boundary and loading conditions of the test were that the inferior surface of L5 vertebra was fixed, and 10 Nm moment and a 150 N pressure preload were applied to the superior surface of L1 vertebra.

Compared with normal model and finest model, the variation of ROM was within 1.03%

in flexion (less than 0.2^o), 4.39% in extension (less than 0.5^o), 0.01% in axial rotation (less than 0.2^o), and 0.001% in lateral bending (less than 0.1^o). From the simulation results, the normal model only required fewer computational times to complete. However, several contact surfaces in facet joint have stress concentration owed to the lower smooth geometry for fewer elements and nodes. Therefore, the appropriate mesh density (finest model) was selected in this study.

Figure 3. 4: Convergence test of the intact model: (A) result of motion changes under flexion;

(B) result of motion changes under extension; (C) result of motion changes under axial rotation;

(D) result of motion changes under lateral bending.

32

3.1.1.2. Follower load setting

During tasks of daily living the lumbar spine withstands compressive loads of very high magnitude along with significant amounts of motion. Compressive loads can easily approach several thousand Newtons during some lifting tasks [104]. Physiologic compressive loads applied to individual lumbar motion segments and the stiffening effect of compressive preload on single functional spine units has been investigated [107]. For the in virtro test, difficulty arises in terms of stability of the lumbar spine when physiologic compressive loads are applied to the entire lumbar spine. The traditional vertical preloads are unable to stabilize the whole lumbar spine specimens under higher physiologic magnitude because the spine without active musculature is unstable at around 100-200 N of vertical preload [108]. The follower load technique, described by Patwardhan et al. [14], applied compressive preload along a path following the lordotic curve of the lumbar spine, and allowed the in vitro spinal models to support higher physiologic loads without damage or instability (Figure 3.5)[14, 109].

Figure 3. 5: (A)The illustration of applying follower load; The experiment setting of follower load ((B) lateral view and(C) front view) [14].

Based on Patwardhan's [14] study, Rohlman [110] used ten fresh-frozen human cadaveric lumbar spines to determine the influence of different loading conditions on

33

intradiscal pressure and intersegmental rotation at all levels of the lumbar spine by applying loads, including pure moments of 3.75, 7.5, and 7.5 Nm plus a follower load of 280 N (140 N on each lateral side)(Figure 3.6). The 280 N corresponds to the partial body weight above the L1 vertebra of a person weighing 66 kg. The results showed that a follower load in the range of the partial body weight is sufficient to stabilize the spine.

Figure 3. 6: Lumbar spine specimen mounted in a spine tester and loaded with a pure moment plus a follower load [110].

Because of the technical limitation of the current experimental set-up, a physiologic compressive preload was applied only while assessing the kinematics in flexion and extension, and was not applied in lateral bending or axial rotation. The preload resulting from muscle activity has a stabilizing effect on a motion segment; therefore, the results pertaining to lateral bending and axial rotation may be viewed as a worst-case scenario.

In finite element studies, the application of a great number of muscle forces is not a problem, but still seldom used as the muscle forces are not known. The influence of muscle forces on the biomechanical behavior of the lumbar spine was investigated using the finite

34

element method. Patwardhan [14] support the following hypothesis concerning the action of muscles in the context of a frontal plane model: Muscle activation causes the internal force resultant to follow a path approximating the tangent of the spinal curve, thereby minimizing the internal shear forces and bending moments and loading the whole lumbar spine in nearly pure compression.

Various numerical models have been used to show the follower load activation effectively on the lumbar spine. Shirazi-Adl [111, 112] developed a method to apply physiologic compression to the lumbar spine through the use of posture changes and “wrapping” elements which wrap around prescribed spatial targets in the center of the endplates of each motion segment such that the compressive load remains perpendicular to the mid-plane of each disc (Fig 3.7). The caudal vertebra S1 is fixed while the remaining vertebrae L1-L5 are unconstrained. Axially fixed compression loads of up to 2800 N are incrementally applied at all L1-L5 vertebral centers (80% at the L1 and the rest evenly distributed among remaining L2-L5 vertebrae to account for differential compression loads along the lumbar spine).

(A) (B)

Figure 3. 7: (A)The illustration of wrapping element; (B)the wrapping element applied on the five levels lumbar spine [111,112].

In Renner's [113] study, A follower load was simulated at each motion segment in the

35

model through a pair of two-node thermo-isotropic truss elements. The follower load trusses were attached bilaterally to the cortical shell of the vertebrae of each motion segment such that each truss spanned the disc, approximately passing through the instantaneous center of rotation of each motion segment, optimizing the follower load path (Figure 3.8). Compressive load was applied to each motion segment by inducing contraction in each of these truss elements by decreasing the temperature in each truss. The results demonstrate that the ability of a large follower load to stiffen the spine in all three planes. Because the follower load concept acts to mimic optimized muscles forces, this study illustrates the important role of muscles in providing spinal stability.

Figure 3. 8: Lateral view of Renner's[113] finite element model showing follower load trusses at each vertebra.

In Rohlmann's [114] study, the weight of the upper body acts in the center of gravity.

Since flexion and extension occur in the sagittal plane, only four muscle groups were simulated: left and right erector spinae, plus left and right rectus abdominis. A compressive follower load was applied to substitute for the unknown stabilizing effect of local muscles

36

(Figure 3.9). The follower load was accomplished by forces of constant magnitude acting in the centers of adjacent vertebral bodies (Figure 3.10).

Applying a follower load instead of a great number of small forces simulating the local dorsal muscles makes realistic loading in in vitro studies feasible. It seems that for sagittal plane motion few global muscle forces are sufficient to achieve realistic results. Rohlmann's [114] results showed that the follower load is a suitable tool to adjust the intradiscal pressure to physiological values without significantly affecting intersegmental motion. However, the local muscles act at a larger lever arm to the center of rotation than a follower load and thus have a greater stabilizing effect. Therefore, local muscles contribute somewhat more to the equilibrium of spinal moments than the follower load.

Figure 3. 9: Finite element model of the lumbar spine with the loads applied in Rolhman's study[114].

37

Figure 3. 10: The illustration of follower load of Rolhman's model [114].

In this study, two kinds of follower load setting for the finite element analysis were proposed. In order to obtain a suitable setting for the follower load. The simpilied follower load was applied on the center of the vertebra (Figure 3.11). In this setting, two hypotheses are used. Firstly, a follower load was applied in the centers of vertebral bodies to substitute for the unknown stabilizing effect of local muscles [114]. Secondly, the thermo-isotropic truss element is used to guarantee the resultant force always toward the center of the next vertebra [113]. The second follower load setting is more close to the experimental set-up. The bilateral follower load was applied through the cable guides on the both side of the spine (Figure 3.12).

The two-node truss element is used to present the cable. The downward force was applied at the end of the cable.

38

Figure 3. 11: The illustration of follower load (simplified: applied in the center).

Figure 3. 12: The illustration of follower load (bilateral: applied bilaterally)

In order to estimate both kinds of follower loads, a 10 Nm moments for four physiological motions were applied. And there are four kinds of loading: 0 N (pure moment), 150 N (vertical load), 400 N (follower load: simplified) and 400 N (follower load: bilateral).

The results are shown in Figure 3.13. The results show that the simplified and the bilateral follower load have similar performance in flexion, extension and rotation. But the bilateral follower load setting shows significantly high stiffness in lateral bending. It is because of the technical limitation of this set-up. Therefore, the results pertaining to lateral bending has the worst correlation to the relistic case. Also, for further research, the position of the lateral plate

39

of the Latero device is going to interfere with the guide device for the cable. Finally, since the bilateral follower load setting is worse, the simipied follower load setting was used in further study.

Figure 3. 13: Range of motion (ROM) calculated for the L1-L5 segments of intact lumbar spine is compared to four kinds of preload condition.

3.1.1.3. Validation of INT model

For the validation of the INT model with follower load, the ROM of the intact model under different loading moments was compared to Rohlmann’s[104] in vitro cadaveric study.

Under 7.5 N-m moments without preload and with a 280 N follower load, the total ROM of five segments lumbar were within one standard deviation in flexion-extension, axial rotation, and lateral bending. The results are shown in Figure 3.14. The present model was verified for further simulations.

40

Figure 3. 14: Range of motion (ROM) calculated for the L1-L5 segments of intact lumbar spine is compared to previous in vitro experiments. (A)Intact lumbar spine without follower load;

(B)intact lumbar spine with simplified follower load.