- Email: [email protected]
A biomechanical investigation of thoracolumbar burst fracture under vertical impact loads using finite element method
Clinical Biomechanics 68 (2019) 29–36
Contents lists available at ScienceDirect
Clinical Biomechanics journal homepage: www.elsevier.com/locate/clinbiomech
A biomechanical investigation of thoracolumbar burst fracture under vertical impact loads using finite element method
T
⁎
Li-Xin Guo , Wu-Jie Li School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Burst fracture Vertical impact load Nonlinear material property Finite element analysis Thoracolumbar spine
Background: A sudden vertical impact load on spine can cause spinal burst fracture, especially in the thoracolumbar junction region. This study aimed at investigating the mechanism of spinal burst fracture under different energy vertical impact loads, producing the failure risk region to understand burst fracture, reducing nervous system damage and guiding clinical treatment. Methods: A nonlinear finite element model of T12-L1 motion segment was created to analyze the response of the vertical impact load. A rigid ball was used to impact the segment vertically to simulate the vertical impact load in practice. There were three different mass balls to represent the different loads: low energy, intermediate energy and high energy (respectively 13 J, 30 J and 56 J). The results of impact force, vertical displacement, stress, intradiscal pressure and contact force were obtained during the process. Findings: At low energy condition, the rigid ball rebounded rapidly. At intermediate energy condition, fractures were initiated in vertebral foramen and left rear regions on the superior cortical bone near the superior endplate of L1. At high energy condition, burst fracture occurred and a part of L1 was isolated from the model. Interpretation: The fracture occurred on the L1 segment only at the intermediate energy and high energy. The strength of vertebral body under low and intermediate energy was enough to support the impact. The burst fracture pattern at high energy was also observed in clinical practice. The findings may explain the mechanism of burst fracture.
1. Introduction The thoracolumbar junction region (TLJ) is a common site of spinal injuries. The fracture of T11-L2 accounts for about half of all spinal fractures (Denis and Perry, 1982; Holdsworth, 1970). Although burst fracture (BF) only comprise about 15% of the TLJ injuries (Denis and Perry, 1982), the incidence of neurological deficit caused by BF is as high as 60% (Denis and Perry, 1982; Mcevoy and Bradford, 1985; Trafton Jr, 1984). So BF of spine is a serious medical condition and impact load is the most frequent cause of trauma. Impact load is often encountered by sportsmen, soldiers and human beings involved extreme activities, like high height falling and traffic accidents. Holdsworth (1970) hypothesized the two-column theory and defined the term “burst fracture” as a fracture that was caused by axial loads. The fracture led to the disc herniation through the upper endplate and resulted in vertebra disruption. Denis and Perry (1982) proposed the concept of middle column from the two-column theory that the anterior and middle column fracture under axial compression loads resulted in the posterior vertebral fragment invading into the spinal
⁎
Corresponding author. E-mail address: [email protected] (L.-X. Guo).
https://doi.org/10.1016/j.clinbiomech.2019.05.018 Received 7 January 2019; Accepted 10 May 2019 0268-0033/ © 2019 Elsevier Ltd. All rights reserved.
canal. Axial compression load is the main reason for BF in fracture classification literature (Magerl et al., 1994). Fractures caused by axial impact loads had been studied in experimental researches (Alpantaki et al., 2010; Daffner and Daffner, 2002; Dai et al., 2007; Fredrickson et al., 1992; Kifune et al., 1995; Kifune et al., 1997; Willén et al., 1984). They focused on the reproduction of BF under different impact energies, the location of BF on the vertebral body, the encroachment of spinal canal fragment, the imaging features of fractured vertebral body and the treatment methods for BF. Cotterill et al. (1987) firstly described a common set-up of drop tower to produce an experimental BF by high energy impact. Germaneau et al. (2017) improved a set-up of Charpy pendulum adapting for an up to 7 vertebrae long segment while imposing pure axial load to reproduce BF. In order to evaluate the severity of BF, some researches (Brandolini et al., 2014; Nicola et al., 2014; Panjabi et al., 1995; Panjabi et al., 1998; Wilcox et al., 2003) measured the changes in the geometry of the spinal canal and parameters of vertebral body, such as, interpedicular widening, vertebral height and vertebral bulge, after axial load impact on the short segment of TLJ. In finite element (FE) studies, Silva et al. (1998) established a FE model of
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
of cancellous bone were with different architecture and mineral density, in this light, the cancellous bone was meshed into seven zones (A–G, Fig. 1) (Zhao and Pollintine, 2009). The intervertebral disc (IVD) was partitioned into two regions, the inner nucleus pulposus and the peripheral annulus fibrosus. The volume of nucleus was set approximately 44%, annulus 56% (White and Panjabi, 1990). The disc was meshed with 8-node hexahedral elements and meshed by 7 longitudinal layers. The homogeneous annulus matrix was reinforced by seven concentric laminar layers of collagen fibers using unidirectional, tension-only 3D spring elements with crosswise pattern at about ± 35° angles in relation to the adjacent endplates (Schmidt et al., 2007). The total volume of fibers was assumed 19% of the annulus ground according defining cross-section of fibers (Natarajan and Andersson, 1999). The six surrounding ligaments, anterior longitudinal ligament (ALL), posterior longitudinal ligament (PLL), ligamentumflavum (LF), capsular ligament (CL), interspinous ligament (ISL) and supraspinous ligament (SSL) (Fig. 1) were modeled in uniform with 1-mm thickness with 3-node (CL only) or 4-node shell elements. The geometric parameters of these ligaments were derived from anatomical literature (Pintar et al., 1992). In total, the mesh of T12-L1 model contained 49,453 nodes and 257,069 elements, ranging from 0.8 mm to 3.0 mm. All elements were set to the first order element. The grid quality was good as shown in supplementary material.
vertebral sections and assigned different materials according to the mineral density from the CT scans to predict the failure loads and fracture patterns. Qiu et al. (2006) built a FE model of T12-L1 to investigate mechanism of vertebral BF depended on the distribution of internal stress. Wilcox et al. (2004) combined experimental results with FE results to study the process of BF. El-Rich et al. (2009) assigned the failure criterion to the vertebral body to reproduce the fractures of L2L3 motion segment under different flexion and extension rates. The aforementioned studies have definitely provided valuable insight to the mechanism of BF. To our best knowledge, few biomechanical studies have been performed to research nonlinear material properties with element rupture under different axial impact loads. The aims of our study are to investigate the mechanism of the burst fracture of TLJ, provide insight into the fracture of TLJ, obtain the variations before and after the fracture and verify the hypothesis in this study that only anterior and middle column fracture during BF in pure axial compression loads. 2. Materials and methods 2.1. The construction of T12-L1 model To establish a three-dimensional (3D) nonlinear complex FE model of T12-L1, 0.6-mm-thick computed tomography (CT) images of spine from a 50th percentile Chinese healthy female(no history of spinal disease) in the lying posture were obtained (Guo and Yin, 2019). The CT images were operated in the software (Mimics 10.0; Materialise Technologies, Leuven) to generate 3D geometry of TLJ. Then the geometry data was imported into CATIA V5R21 (Dassault Systèmes, VélizyVillacoublay) to create the 3D solid model of T12-L1. The model was meshed in pre-processing software of ANSA (BETA CAE Systems). The cortical bone and endplates modeled by 3-node shell elements consisted of the outer vertebral body and outer posterior elements. The interparts then were filled with cancellous bone core meshed by 4-node solid element. According to local variations of cortical and endplate thickness, the outer bone was subdivided into nine different thickness zones (Fazzalari et al., 2006; Roberts et al., 1997) (Fig. 1). Various locations
2.2. Material property All bony components were assumed to be homogenous and followed a strain rate-dependent elastoplastic Johnson-Cook material law (Table 1) that allowed computing von-Mises hardening with ductile damage until potential rupture (El-Rich et al., 2009). Before the equivalent stress reached to the yield stress, the material behaved as linear elastic. For higher stress values, the material behavior was plastic and the equivalent stress at constant temperature was calculated using the following Eq. (1):
σ = (A + Bε n)(1 + Clnε ∗̇ )
(1)
Fig. 1. The finite element model ofT12-L1 and Subdivision of the cortical bone, the cancellous bone into 9 components with different cortical thickness and into 7 components with different material properties. 30
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Table 1 Summary of the material properties of bone and IVD tissues used in the model (adapted from the reference (El-Rich et al., 2009; Lee et al., 2000; Schmidt et al., 2007; Wagnac et al., 2012)). Cancellous (A–G zone) and endplate centera
Material properties
Density (×10−6 kg/mm3) Modulus of elasticity, E (MPa) Poisson ratio, ν Yield stress, A (MPa) Harding modulus, B (MPa) Hardening exponent, n Failure plastic strain, εp Maximum stress (MPa) Strain rate coefficient, C Reference strain rate, ε0̇
Material properties
−6
A
B
C
D
E
F
G
1.8 93.7 0.25 1.95 8.5 1 0.082 2.65 0.533 0.008
1.8 93.7 0.25 1.95 7.0 1 0.06 2.3 0.533 0.008
2.0 93.7 0.25 1.95 8.5 1 0.082 2.65 0.533 0.008
2.0 93.7 0.25 1.95 8.1 1 0.08 2.6 0.533 0.008
2.5 93.7 0.25 1.95 12.5 1 0.104 3.25 0.533 0.008
2.5 93.7 0.25 1.95 12.5 1 0.104 3.25 0.533 0.008
1.8 93.7 0.25 1.95 7.0 1 0.06 2.3 0.533 0.008
Nucleus pulposus
3
Density (×10 kg/mm ) Poisson ratio, ν C10(MPa) C01(MPa) a
Cortical and endplate margin
0.2 4014 0.3 142.5 1338 1 0.071 190 0.272 0.008
Annulus groundsubstance
Low strain rate
High strain rate
Low strain rate
High strain rate
Fibers
1 0.495 0.12 0.03
1 0.495 31.8 8.0
1.2 0.45 0.18 0.045
1.2 0.45 11.8 2.9
Nonlinear
Bony endplate center material properties are the same as zone B.
Fig. 2. (A) Loading and boundary conditions of impact, (B) Loading and boundary conditions of quasi-static compressive test of the disc, (C) Compressive forcedisplacement curve of the disc, (D) loading and boundary conditions for validation of moment in a standard 3-D coordinate system, (E) degrees of left and right axial rotation, (F) degrees of flexion and extension, (G) Degrees of left and right lateral bending.
ε /̇ ε ̇ ε ̇ ≥ ε0̇ ε ∗̇ = ⎧ 0 ⎨ ⎩ 1 ε ̇ < ε0̇
Here, C10, C01 material constants; I1 and I2 the first and second invariants of the deviatoric component of the left Cauchy-Green strain tensor; J = V/V0, local volume ratio; d = 2/K, material incompressibility factor, K initial bulk modulus. But the law does not consider strain rate which has a significant influence on the mechanical behavior of IVD. Two sets of material parameters were adopted according to different strain rates (Schmidt et al., 2007; Wagnac et al., 2011). One set named as low strain rate was used in static or quasi static condition and the other set named high strain rate was used in high strain rate dynamic condition, such as the impact in this study. The property of collagen fibers was nonlinear based on a stressstrain curve (Shirazi-Adl et al., 1986). Since the stiffness and crosssection of fibers increase from the inner layers to the outer layers of the annulus, the revised ratio of the cross-sectional areas and elasticity
(2)
where, σ equivalent stress, A yield stress, B hardening modulus, n hardening exponent, ε ̇ current strain rate, ε0̇ reference strain rate and εpplastic strain (true strain). When εpreachesεmax, either for tension, compression or shear, shell elements would be deleted and the solid elements are not deleted but the deviatoric stress of solid elements are permanently set to zero. The IVD was governed by the hyper-elastic material law (Table 1) based on a first-order Mooney-Rivlin formulation with strain energy density equation, W is defined as follows.
W = C10 (I1 − 3) + C01 (I2 − 3) + (J − 1)/ d
(3)
K = 6(C10 + C01)/3(1 − 2ν)
(4) 31
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Fig. 3. The simulation results. (A)impact forces between the rigid ball and the rigid plate under different energy conditions, (B)the vertical displacements of centroid of the superior surface of T12 under different energy conditions, (C) Total contact forces between the facet surfaces under different energy conditions, (D) IDP changes under different energy conditions.
possible rebound (not impact tie contact). The ball was constrained and allowed only vertical motion to ensure vertical impact for minimizing the undesired flexion or extension. The inferior surface and spinous process of L1 were fixed in all degrees. The Radioss block100 explicit dynamic solver (Altair Hyperworks Inc. USA) was used to solve the investigation. The solution time was set to 6 ms under HE condition and 8 ms under LE and IME conditions and the running duration was about 4 h (Computer Configuration: Intel i52400, 3.10 Hz, 4GB RAM, Windows 7).
constants of the seven layers were adopted (Shirazi-Adl et al., 1986). The properties of ligaments were modeled in nonlinear material properties (Sharma et al., 1995). 2.3. Contact, loading and boundary conditions Due to the different mesh methods of bone, IVD and ligaments, a tied contact was created among them. Frictionless contact interfaces were used between adjacent facet joints to compute the contact forces. Frictionless contact interfaces were also created between different parts to avoid any possible penetration. All the loading and boundary conditions were set according to the literature (Kifune et al., 1995; Kifune et al., 1997; Panjabi et al., 1995; Panjabi et al., 1998; Qiu et al., 2006) (Fig. 2A). According to their experimental results, the TLJ initiated fracture in the impact energy of 2 kg ball (30 J) and BF in 4.3 kg (56 J). These were corresponded to the results that the energies replicating BF with no modification of the spinal dynamics are around 50 and 180 J (Germaneau et al., 2017; Jones et al., 2011; Panjabi et al., 2000). 1 kg, 2.3 kg and 4.3 kg were assigned to the rigid shell ball to represent low energy (LE), intermediate energy (IME) and high energy (HE) impact. In order to transfer the impact loads down to T12 uniformly and synchronously, a 3-mm thick rigid shell plate was created by the outline of T12 top surface. The shell plate was tied to the upper surface of T12. The rigid shell ball with the radius of 11 mm was placed over the rigid plate following its center vertically above the centroid of the plate with a distance of 5 mm between the shell plate and the lowest point of the ball. The ball was set at a 5.1 m/s downward velocity without gravitational acceleration to replace the experimental velocity of free falling from the height of 1.4 m for reducing computational cost. Contact was created between the rigid ball and rigid plate with
2.4. Validation of the model A quasi-static load was applied to the model of TLJ excluding posterior elements according to previous references (El-Rich et al., 2009; Zeng et al., 2013) (Fig. 2B). The loading rate was 1.26 mm/s and the loading duration was 1.2 s. KEREL order (kinetic energy relaxation) had been used to minimize the effect of kinetic energy and a curve of loaddisplacement was obtained. The curve was compared with the reported results (El-Rich et al., 2009; Zeng et al., 2013). The present model was validated under all physiological loading modes- flexion/extension, left/right lateral bending, and left/right axial rotation, by comparing these results with experimental results (Markolf, 1972; Oxland et al., 1992) and numerical results (Qiu et al., 2006) according to the range of motion (RoM) of the TLJ. The 7.5 Nm pure moment was applied to the TLJ in 3 incremental steps, 2.5 Nm, 5 Nm and 7.5 Nm, to create the relevant motion (Fig. 2D). 3. Results A compressive force-displacement curve of the 3-D nonlinear FE model of T12-L1was obtained by a velocity loading as showed in 32
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
initiated in vertebral foramen on the left-rear region on the superior cortical bone near the superior endplate of L1 at 1.85 ms. Then the fractures propagated along the circumferential direction. At 2.05 ms and 2.20 ms, fractures also occurred in the right-rear of cortical bone near the superior endplate and the posterior region of superior endplate of L1. The stress on cancellous bone adjacent to the fractural cortical bone was found to be zero meaning that the fracture occurred in the region. The fractural time was a little early than the cortical bone at 1.55 ms. Several fractural cancellous bone elements were located in the middle region of L1 on the left at the same time. The final fractural cancellous bone represented by dark blue elements was shown in Fig. 4C. In HE condition, more severe fracture occurred. The locations and the time of initial fractures were the same as those in IME, but the range of fracture was larger and the following propagation speed faster. At 2.8 ms, the fracture occurred in the middle of the antero-superior part of L1. With the rigid ball continuously compressing, the fractures in different locations connected together and formed a circle which divided L1 column into two parts. The stress distribution of cancellous bone was similar to that of IME, but the fracture was more severe. The superior and the anterior part of L1 collapsed toward the core of cancellous bone (Fig. 4D). As shown in Figs. 4 and 5, the cortical elements in the high stress region (red zone) exceeding the failure plastic strain would be deleted suddenly. The color of cancellous elements in the high stress region exceeding the failure plastic strain would turn to dark blue, but the elements keep deformation with zero stress output.
Fig. 2C. The curve fell within the quasi-static load-displacement experimental and FE corridors and showed a nonlinear compressive behavior. The predicted load-angular curves under various pure moments showed similar trends of incremental rotational motion and matched well with the published results in Fig. 2E–G. The validation study demonstrated that the present model was realistic. The impact forces between the rigid plate and the ball were gained (Fig. 3A). The HE reached to the maximum impact force of about 13KN. The minimum impact force was about 10.5 KN for LE. The vertical displacement was recorded by the centroid of the superior surface of T12 (Fig. 3B). According to Fig. 3A, the ball of LE first rebounded at about 1.8 ms and then the ball of IME rebounded at about 3.5 ms. The HE ball maintained a vertical motion without rebound. The maximum intradiscal pressure (IDP) caused by the three impacting energies were 44.5 MPa,50.2 MPa and 51.0 MPa respectively under the maximum impact forces (Fig. 3C). The rebound caused the IDP decreasing and vibrating in LE. The same case also happened in IME. The fractures in the cortical bone and endplate also decreased the IDP. The negative values of IDP represented the state of tension under impacting loads. The response from total contact force was markedly different under three impacting energies (Fig. 3D). In LE, the facet contact force exhibited a periodic motion, which was a cyclic contact between adjacent facets. The contact forces of IME and HE showed a vibrating attenuate tendency after the maximum impact force. The adjacent facets kept contact during the HE impacting time. The adjacent facets separated after the rebound and started a cyclic contact (Fig. 3D) under LE condition. In LE condition, the LE did not cause fracture of vertebral body. High stresses within cortical bone were located in the left and right rear region of L1 near the superior endplate margin and between the pedicles. The maximum stress reached to 237 MPa (Fig. 4A). High stress on cancellous bone was in the same region and the maximum stress was 11 MPa. In IME condition, the burst fracture occurred. Fractures were
4. Discussion This model used an accurate and practical method to investigate the mechanism of BF including the initiation, propagation and termination. This study not only defined high risk region owing to high-stress concentration (Lee et al., 2000; Qiu et al., 2006; Wilcox et al., 2004), but also reproduced the process of burst fracture. The load of vertical
Fig. 4. (A) The von Mises stress distribution in cortical bone under LE condition corresponding to the maximum vertical displacement of T12 (excluding the posterior part), (B) Initiation and propagation of cortical fracture in L1 under IME condition, (C) Initiation and termination of fracture of cancellous bone in L1 under IME condition, (D) Initiation and propagation of cortical fracture in L1 under HE condition. 33
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Fig. 5. The stress distribution of burst fracture of L1 cancellous bone under HE condition. (A) to (D) The fractures of cancellous bone at different time in different views, The time was set at 1.8 ms, 2.0 ms, 2.1 ms, 2.3 ms, 3.3 ms and 5 ms. (A) the left view, (B) the right view, (C) the front view, (D) the sagittal plane.
Depending on the failure strength of these regions, failure could be initiated in any of them. In this study, the initiation of fractures was on the posterior edge of superior endplate of L1 and adjacent cortical bone (including posterior cortical wall) (Fig. 3C and Fig. 4). The three results had shown varying degrees of damage. In the LE impacting condition, the TLJ was intact and no fracture was noted on the vertebral body. In the IME impacting (Fig. 4B–C), fractures occurred in both the cancellous bone and the cortical bone. The fractures initiated on the left posterior edge of superior endplate of L1 and the posterosuperior wall of vertebral body between the pedicle bases. Then the initial fractures trended to connect together and there was a new fracture on the right posterior edge of the superior endplate of L1. The fractures of cancellous bone adjacent to the fractured cortical bone occurred in the same time. A fractural zone in the middle of anterior column was observed where a potential fractured zone of burst fracture was described in the previous literature (Denis and Perry, 1982; Hakim and King, 1979; Panjabi et al., 1995). The fractures spread toward the core of the cancellous bone internally and would form an oblique line in the sagittal plane. Although the fractures occurred in IME, the vertebral body was a whole entity which had enough stiffness to rebound the impacting ball. The fragments of posterior wall invading the spinal cord recovered as before. There was no obvious fracture-dislocation of the TLJ. The structure of TLJ recovered after the impacting ball left the vertebral body. In HE (Fig. 5), the first half of fracture was similar to the fracture under the IME impacting condition. In the second half, with the impacting ball continuously compressing, the fracture continued to expand. Compared with the result of IME, a part of the anterior cortical bone and a part of the superior cancellous bone were isolated from the vertebral body along the fractural line in the sagittal plane resulting the
impacting ball was a transient variation with high to low velocity, which was closer to the fact in daily life but was different to the constant load (El-Rich et al., 2009; Fradet et al., 2014; Mustafy et al., 2014; Wagnac et al., 2011). The posterior ligament complex is seldom ruptured under axial compression unless the spine is subjected to forcible rotation or flexion-rotation (Holdsworth, 1970). The failure model of posterior ligament complex was simplified and was not analyzed in this study. Many researchers have offered classification and extensive description of fracture patterns in clinical observation. Magerl et al. (1994) described a total of 27 fracture patterns based on different single loads or combining forces. Vaccaro et al. (2005) evaluated the severity of spinal injury by measuring the parameters of spine. The FE method was a useful tool to verify the fracture patterns and provide the spinal injury severity score. Roaf (1960), Tran et al. (1995) and Qiu et al. (2006) suggested that increased pressure in the nucleus pulposus caused endplates to bulge toward the cancellous core resulting to crack which allowed the nucleus material to enter the vertebral body at a fast rate than the content could be expelled. The three results of this study appeared to contradict their hypotheses. The IDP only increased rapidly with the increased impact force. As the impact force decreased, the IDP in the nucleus never reached a maximum peak any longer. The bulge of IVD in the transverse direction was larger than it in the vertical direction meaning the rate of entering into the endplate was not great. This was consistent with the finding of Wilcox et al. (2004). The pressure and stress in the center of IVD were lower than those on the periphery which also supported the opinion. Lee et al. (2000) reported that high stresses could lead to fracture in the endplate and the posterior surface of the cortical wall.
Fig. 6. The patterns of burst fracture. (A)The FE result of this study, (B)Type B of burst fracture classified by Denis et al., (C)Type A of burst fracture classified by Magerl et al., (D)Lateral tomogram of a burst fracture. 34
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
and after burst fracture on the special model. The hypothesis that only anterior and middle columns fracture during BF in pure axial compression load was verified in this study. The burst fracture type under the HE condition was observed in clinical practice. The study suggested an integrated, detailed vision and insight into the mechanism, instability and treatment of TLJ burst fracture in clinical and experimental studies.
spinal cord stenosis. The vertebral body of L1 was divided into two parts (Fig. 6A). The same burst fracture was classified as type B by Denis and Perry (1982), which was described as disruption of superior endplate, loss of height, fracture of posterior wall of vertebral body (Fig. 6B and D)and was classified as type A3.1 by Magerl et al. (1994), which the upper or lower half of the vertebral body has burst while the other half remains intact (Fig. 6C). If the impact is not strong, the bone is within the elastic deformation, in ideal circumstances, that the bone will absorb the energy and start the reciprocating vibration, the ball can rebound. If the impact is strong, the plastic deformation will occur. The bone will absorb the energy and not release the energy which caused the plastic deformation. In this study, the ball rebound under LE and IME conditions. In LE condition, the velocity of rebound was 4.03 m/s, 5 J was absorbed by model. In IME condition, the velocity of rebound was 2.2 m/s, 9.8 J was absorbed by model. In HE condition, the plastic deformation was too large to converge in the analysis. Facet joints could play an important role during the impacting process. A hypothesis proposed by Wilcox et al. (2004) that the facet joints transmitted shear forces to the posterosuperior part of vertebral body. The fact that higher stress around the pedicle bases with high contact forces supported this hypothesis. The maximum IDP values in IME and HE were approximate which was different in LE. This may attribute to the sudden penetration of the nucleus into the gap of fractural bone, releasing the pressure. In this study, the trajectory of the ball is orthogonal to the plane of the endplate, so the rigid ball would not produce unexpected moment. The additional moment change the facet joint force and the morphology of IVD which have the contribution of different BF. This was discussed in Langrana et al. (2002) and had been completed in our following study. The numerical simulation has its limitations. It has been reported that the unique geometrical features characterizing each spine can determine significant changes in load distributions within the spine segments, so the lack of the diversity in resulting fracture patterns, various fractures under vertical impacting classified by the researchers (Denis and Perry, 1982; Holdsworth, 1970; Magerl et al., 1994) could not be obtained. This study could not compare between all the compressive fracture types. The study did not model the active and passive muscle tissues which may have influence on the TLJ kinematic responses. In high velocity dynamic loading, it was reasonable that there was no fluid flow in spinal segments as the aim of this study was focused on the mechanism of burst fracture under vertical impacting. It was different from the long-term (Argoubi and Shirazi-Adl, 1996; Schroeder et al., 2010; Heuer et al., 2007) compressive loads which needed visco-poroelastic property to study in a time frame of minutes or hours. But it still could be improved in a limited fluid-flow material property and independent mechanism under high velocities (Disilvestro and Suh, 2001; Oloyede and Broom, 1993). The load in this study was only axial compression and this could be improved by applying combined or complex loads to obtain more fracture patterns in clinical observations. The validation guaranteed the intended usage of FE model. The use of a verified model in static bending tests, especially, in quasi-static tests is one of the strengths of this study, but it still could be improved in dynamic validation which had a diverse response compered to in static condition. The L1-L5 dynamic experimental test results were obtained by Stokes and Gardner-Morse (2016), the TLJ segment experimental data was not found, and it should be completed in future work.
Declaration of Competing Interest The authors declare that they have no conflict of interest. Acknowledgement This work was supported by the National Natural Science Foundation of China (51875096, 51275082). References Alpantaki, K., Bano, A., Pasku, D., Mavrogenis, A.F., Papagelopoulos, P.J., Sapkas, G.S., Korres, D.S., Katonis, P., 2010. Thoracolumbar burst fractures: a systematic review of management. Orthopedics 33, 422–429. Argoubi, M., Shirazi-Adl, A., 1996. Poroelastic creep response analysis of a lumbar motion segment in compression. J. Biomech. 29, 1331–1339. Brandolini, N., Kapur, N., Hall, R.M., 2014. Dynamics of interpedicular widening in spinal burst fractures: an in vitro investigation. Spine J. 14, 2164–2171. Cotterill, P.C., Kostuik, J.P., Wilson, J.A., Fernie, G.R., Maki, B.E., 1987. Production of a reproducible spinal burst fracture for use in biomechanical testing. J. Orthop. Res. 5, 462–465. Daffner, R.H., Daffner, S.D., 2002. Vertebral injuries: detection and implications. Eur. J. Radiol. 42, 100–116. Dai, L.Y., Jiang, S.D., Wang, X.Y., Jiang, L.S., 2007. A review of the management of thoracolumbar burst fractures. Surg. Neurol. 67, 221–231. Denis, Francis, Perry, John F., 1982. The three-column spine concept and its significance in the classification of thoracolumbar spinal injuries. J. Trauma 22. Disilvestro, M.R., Suh, J.K., 2001. A cross-validation of the biphasic poroviscoelastic model of articular cartilage in unconfined compression, indentation, and confined compression. J. Biomech. 34, 519–525. El-Rich, M., Arnoux, P.J., Wagnac, E., Brunet, C., Aubin, C.E., 2009. Finite element investigation of the loading rate effect on the spinal load-sharing changes under impact conditions. J. Biomech. 42, 1252–1262. Fazzalari, N.L., Parkinson, I.H., Fogg, Q.A., Sutton-Smith, P., 2006. Antero-postero differences in cortical thickness and cortical porosity of T12 to L5 vertebral bodies. Joint Bone Spine 73, 293–297. Fradet, L., Petit, Y., Wagnac, E., Aubin, C.E., Arnoux, P.J., 2014. Biomechanics of thoracolumbar junction vertebral fractures from various kinematic conditions. Med. Biol. Eng. Comput. 52, 87–94. Fredrickson, B.E., Edwards, W.T., Rauschning, W., Bayley, J.C., Yuan, H.A., 1992. Vertebral burst fractures: an experimental, morphologic, and radiographic study. Spine 17, 1012–1021. Germaneau, A., Vendeuvre, T., Saget, M., Doumalin, P., Dupré, J.C., Brémand, F., Hesser, F., Brèque, C., Maxy, P., Roulaud, M., 2017. Development of an experimental model of burst fracture with damage characterization of the vertebral bodies under dynamic conditions. Clin. Biomech. 49, 139–144. Guo, L.X., Yin, J.Y., 2019. Finite element analysis and design of an interspinous device using topology optimization. Med. Biol. Eng. Comput. 57, 89–98. Hakim, N.S., King, A.I., 1979. A three dimensional finite element dynamic response analysis of a vertebra with experimental verification. J. Biomech. 12, 277–292. Heuer, F., Schmitt, H., Schmidt, H., Claes, L., Wilke, H.J., 2007. Creep associated changes in intervertebral disc bulging obtained with a laser scanning device. Clin. Biomech. 22, 737–744. Holdsworth, F.W., 1970. HoldsworthFFractures, dislocations, and fracture-dislocations of the spine. J Bone Joint Surg 52A: 1534-1551. J. Bone Joint Surg. 52, 1534–1551. Jones, H.L., Crawley, A.L., Noble, P.C., Schoenfeld, A.J., Weiner, B.K., 2011. A novel method for the reproducible production of thoracolumbar burst fractures in human cadaveric specimens. Spine J. 11, 447–451. Kifune, M., Panjabi, M.M., Arand, M., Liu, W., 1995. Fracture pattern and instability of thoracolumbar injuries. Eur. Spine J. 4, 98–103. Kifune, M., Panjabi, M.M., Liu, W., Arand, M., Vasavada, A., Oxland, T., 1997. Functional morphology of the spinal canal after endplate, wedge, and burst fractures. J. Spinal Disord. 10, 457–466. Langrana, N.A., Harten, R.D.R.D., Lin, D.C., Reiter, M.F., Lee, C.K., 2002. Acute thoracolumbar burst fractures: a new view of loading mechanisms. Spine 27, 498–508. Lee, C.K., Kim, Y.E., Lee, C.S., Hong, Y.M., Jung, J.M., Goel, V.K., 2000. Impact response of the intervertebral disc in a finite-element model. Spine 25, 2431–2439. Magerl, F., Aebi, M., Gertzbein, S.D., Harms, J., Nazarian, S., 1994. A comprehensive classification of thoracic and lumbar injuries. Eur. Spine J. 3, 184–201. Markolf, K.L., 1972. Deformation of the thoracolumbar intervertebral joints in response to external loads: a biomechanical study using autopsy material. J. Bone Joint Surg. (Am. Vol.) 54, 511–533.
5. Conclusion A detailed FE model of T12-L1 with nonlinear material property and failure criterion was built to research the mechanism of burst fracture under vertical impacting loads. The results provided the perspective on BF including the initiation, propagation and termination in detail, and the variations of stress, IDP, contact force and vertebral body before 35
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Pa 1976) 11, 914–927. Silva, M.J., Keaveny, T.M., Hayes, W.C., 1998. Computed tomography-based finite element analysis predicts failure loads and fracture patterns for vertebral sections. J. Orthop. Res. 16, 300–308. Stokes, I.A.F., Gardner-Morse, M., 2016. A database of lumbar spinal mechanical behavior for validation of spinal analytical models. J. Biomech. 49, 780–785. Trafton Jr., P.G., 1984. BC: computed tomography of thoracic and lumbar spine injuries. J. Trauma 24, 506–515. Tran, N.T., Watson, N.A., Tencer, A.F., Ching, R.P., Anderson, P.A., 1995. Mechanism of the burst fracture in the thoracolumbar spine. The effect of loading rate. Spine 20, 1984–1988. Vaccaro AR, Lehman RA, R John H, Anderson PA, Mitchel H, Rune H, James H, Marcel D, Kirkham W, Fehlings MG: A new classification of thoracolumbar injuries: the importance of injury morphology, the integrity of the posterior ligamentous complex, and neurologic status. Spine 2005, 30:2325. Wagnac, E., Arnoux, P.J., Garo, A., Elrich, M., Aubin, C.E., 2011. Calibration of hyperelastic material properties of the human lumbar intervertebral disc under fast dynamic compressive loads. J. Biomech. Eng. 133, 419–426. Wagnac, E., Arnoux, P.J., Garo, A., Aubin, C.E., 2012. Finite element analysis of the influence of loading rate on a model of the full lumbar spine under dynamic loading conditions. Med. Biol. Eng. Comput. 50, 903–915. White, A.A., Panjabi, M.M., 1990. Clinical Biomechanics of the Spine, 2nd edn. Lippincott, Philadelphia. Wilcox, R.K., Boerger, T.O., Allen, D.J., Barton, D.C., Limb, D., Dickson, R.A., Hall, R.M., 2003. A dynamic study of thoracolumbar burst fractures. J. Bone Joint Surg. Am. 85, 2184–2189. Wilcox, R.K., Allen, D.J., Hall, R.M., Limb, D., Barton, D.C., Dickson, R.A., 2004. A dynamic investigation of the burst fracture process using a combined experimental and finite element approach. Eur. Spine J. 13, 481–488. Willén, J., Lindahl, S., Irstam, L., Aldman, B., Nordwall, A., 1984. The thoracolumbar crush fracture. An experimental study on instant axial dynamic loading: the resulting fracture type and its stability. Spine 9, 624. Zeng, Z.L., Zhu, R., Li, S.Z., Yu, Y., Wang, J.J., Jia, Y.W., Chen, B., Cheng, L.M., 2013. Formative mechanism of intracanal fracture fragments in thoracolumbar burst fractures: a finite element study. Chin. Med. J. 126, 2852–2858. Zhao, F., Pollintine, P., 2009. Bd, Adams M, Dolan P: vertebral fractures usually affect the cranial endplate because it is thinner and supported by less-dense trabecular bone. Bone 44, 372–379.
Mcevoy, R.D., Bradford, D.S., 1985. The management of burst fractures of the thoracic and lumbar spine. Experience in 53 patients. Spine 10, 631–637. Mustafy, T., Elrich, M., Mesfar, W., Moglo, K., 2014. Investigation of impact loading rate effects on the ligamentous cervical spinal load-partitioning using finite element model of functional spinal unit C2-C3. J. Biomech. 47, 2891–2903. Natarajan, R.N., Andersson, G.B., 1999. The influence of lumbar disc height and crosssectional area on the mechanical response of the disc to physiologic loading. Spine 24, 1873–1881. Nicola, B., Nikil, K., Hall, R.M., 2014. Dynamics of interpedicular widening in spinal burst fractures: an in vitro investigation. Spine J. 14, 2164–2171. Oloyede, A., Broom, N.D., 1993. A physical model for the time-dependent deformation of articular cartilage. Connect. Tissue Res. 29, 251–261. Oxland, T.R., Lin, R.M., Panjabi, M.M., 1992. Three-dimensional mechanical properties of the thoracolumbar junction. J. Orthop. Res. 10, 573–580. Panjabi, M.M., Kifune, M., Wen, L., Arand, M., Oxland, T.R., Lin, R.M., Yoon, W.S., Vasavada, A., 1995. Dynamic canal encroachment during thoracolumbar burst fractures. J. Spinal Disord. 8, 39–48. Panjabi, M.M., Kifune, M., Liu, W., Arand, M., Vasavada, A., Oxland, T.R., 1998. Graded thoracolumbar spinal injuries: development of multidirectional instability. Eur. Spine J. 7, 332–339. Panjabi, M.M., Hoffman, H., Kato, Y., Cholewicki, J., 2000. Superiority of incremental trauma approach in experimental burst fracture studies. Clin. Biomech. 15, 73–78. Pintar, F.A., Yoganandan, N., Myers, T., Elhagediab Jr., A., 1992. SA: biomechanical properties of human lumbar spine ligaments. J. Biomech. 25, 1351–1356. Qiu, T.X., Tan, K.W., Lee, V.S., Teo, E.C., 2006. Investigation of thoracolumbar T12-L1 burst fracture mechanism using finite element method. Med. Eng. Phys. 28, 656–664. Roaf, R., 1960. A study of the mechanics of spinal injuries. Bone Joint J. 42-B, 810–823. Roberts, S., Mccall, I.W., Menage, J., Haddaway, M.J., Eisenstein, S.M., 1997. Does the thickness of the vertebral subchondral bone reflect the composition of the intervertebral disc? Eur. Spine J. 6, 385–389. Schmidt, H., Kettler, A., Heuer, F., Simon, U., Claes, L., Wilke, H.J., 2007. Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading. Spine 32, 748–755. Schroeder, Y., Huyghe, J.M., Donkelaar, C.C.V., Ito, K., 2010. A biochemical/biophysical 3D FE intervertebral disc model. Biomech. Model. Mechanobiol. 9, 641–650. Sharma, M., Langrana, N.A., Rodriguez, J., 1995. Role of ligaments and facets in lumbar spinal stability. Spine 20, 887–900. Shirazi-Adl, A., Ahmed, A.M., Shrivastava, S.C., 1986. Mechanical response of a lumbar motion segment in axial torque alone and combined with compression. Spine (Phila
36
Contents lists available at ScienceDirect
Clinical Biomechanics journal homepage: www.elsevier.com/locate/clinbiomech
A biomechanical investigation of thoracolumbar burst fracture under vertical impact loads using finite element method
T
⁎
Li-Xin Guo , Wu-Jie Li School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Burst fracture Vertical impact load Nonlinear material property Finite element analysis Thoracolumbar spine
Background: A sudden vertical impact load on spine can cause spinal burst fracture, especially in the thoracolumbar junction region. This study aimed at investigating the mechanism of spinal burst fracture under different energy vertical impact loads, producing the failure risk region to understand burst fracture, reducing nervous system damage and guiding clinical treatment. Methods: A nonlinear finite element model of T12-L1 motion segment was created to analyze the response of the vertical impact load. A rigid ball was used to impact the segment vertically to simulate the vertical impact load in practice. There were three different mass balls to represent the different loads: low energy, intermediate energy and high energy (respectively 13 J, 30 J and 56 J). The results of impact force, vertical displacement, stress, intradiscal pressure and contact force were obtained during the process. Findings: At low energy condition, the rigid ball rebounded rapidly. At intermediate energy condition, fractures were initiated in vertebral foramen and left rear regions on the superior cortical bone near the superior endplate of L1. At high energy condition, burst fracture occurred and a part of L1 was isolated from the model. Interpretation: The fracture occurred on the L1 segment only at the intermediate energy and high energy. The strength of vertebral body under low and intermediate energy was enough to support the impact. The burst fracture pattern at high energy was also observed in clinical practice. The findings may explain the mechanism of burst fracture.
1. Introduction The thoracolumbar junction region (TLJ) is a common site of spinal injuries. The fracture of T11-L2 accounts for about half of all spinal fractures (Denis and Perry, 1982; Holdsworth, 1970). Although burst fracture (BF) only comprise about 15% of the TLJ injuries (Denis and Perry, 1982), the incidence of neurological deficit caused by BF is as high as 60% (Denis and Perry, 1982; Mcevoy and Bradford, 1985; Trafton Jr, 1984). So BF of spine is a serious medical condition and impact load is the most frequent cause of trauma. Impact load is often encountered by sportsmen, soldiers and human beings involved extreme activities, like high height falling and traffic accidents. Holdsworth (1970) hypothesized the two-column theory and defined the term “burst fracture” as a fracture that was caused by axial loads. The fracture led to the disc herniation through the upper endplate and resulted in vertebra disruption. Denis and Perry (1982) proposed the concept of middle column from the two-column theory that the anterior and middle column fracture under axial compression loads resulted in the posterior vertebral fragment invading into the spinal
⁎
Corresponding author. E-mail address: [email protected] (L.-X. Guo).
https://doi.org/10.1016/j.clinbiomech.2019.05.018 Received 7 January 2019; Accepted 10 May 2019 0268-0033/ © 2019 Elsevier Ltd. All rights reserved.
canal. Axial compression load is the main reason for BF in fracture classification literature (Magerl et al., 1994). Fractures caused by axial impact loads had been studied in experimental researches (Alpantaki et al., 2010; Daffner and Daffner, 2002; Dai et al., 2007; Fredrickson et al., 1992; Kifune et al., 1995; Kifune et al., 1997; Willén et al., 1984). They focused on the reproduction of BF under different impact energies, the location of BF on the vertebral body, the encroachment of spinal canal fragment, the imaging features of fractured vertebral body and the treatment methods for BF. Cotterill et al. (1987) firstly described a common set-up of drop tower to produce an experimental BF by high energy impact. Germaneau et al. (2017) improved a set-up of Charpy pendulum adapting for an up to 7 vertebrae long segment while imposing pure axial load to reproduce BF. In order to evaluate the severity of BF, some researches (Brandolini et al., 2014; Nicola et al., 2014; Panjabi et al., 1995; Panjabi et al., 1998; Wilcox et al., 2003) measured the changes in the geometry of the spinal canal and parameters of vertebral body, such as, interpedicular widening, vertebral height and vertebral bulge, after axial load impact on the short segment of TLJ. In finite element (FE) studies, Silva et al. (1998) established a FE model of
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
of cancellous bone were with different architecture and mineral density, in this light, the cancellous bone was meshed into seven zones (A–G, Fig. 1) (Zhao and Pollintine, 2009). The intervertebral disc (IVD) was partitioned into two regions, the inner nucleus pulposus and the peripheral annulus fibrosus. The volume of nucleus was set approximately 44%, annulus 56% (White and Panjabi, 1990). The disc was meshed with 8-node hexahedral elements and meshed by 7 longitudinal layers. The homogeneous annulus matrix was reinforced by seven concentric laminar layers of collagen fibers using unidirectional, tension-only 3D spring elements with crosswise pattern at about ± 35° angles in relation to the adjacent endplates (Schmidt et al., 2007). The total volume of fibers was assumed 19% of the annulus ground according defining cross-section of fibers (Natarajan and Andersson, 1999). The six surrounding ligaments, anterior longitudinal ligament (ALL), posterior longitudinal ligament (PLL), ligamentumflavum (LF), capsular ligament (CL), interspinous ligament (ISL) and supraspinous ligament (SSL) (Fig. 1) were modeled in uniform with 1-mm thickness with 3-node (CL only) or 4-node shell elements. The geometric parameters of these ligaments were derived from anatomical literature (Pintar et al., 1992). In total, the mesh of T12-L1 model contained 49,453 nodes and 257,069 elements, ranging from 0.8 mm to 3.0 mm. All elements were set to the first order element. The grid quality was good as shown in supplementary material.
vertebral sections and assigned different materials according to the mineral density from the CT scans to predict the failure loads and fracture patterns. Qiu et al. (2006) built a FE model of T12-L1 to investigate mechanism of vertebral BF depended on the distribution of internal stress. Wilcox et al. (2004) combined experimental results with FE results to study the process of BF. El-Rich et al. (2009) assigned the failure criterion to the vertebral body to reproduce the fractures of L2L3 motion segment under different flexion and extension rates. The aforementioned studies have definitely provided valuable insight to the mechanism of BF. To our best knowledge, few biomechanical studies have been performed to research nonlinear material properties with element rupture under different axial impact loads. The aims of our study are to investigate the mechanism of the burst fracture of TLJ, provide insight into the fracture of TLJ, obtain the variations before and after the fracture and verify the hypothesis in this study that only anterior and middle column fracture during BF in pure axial compression loads. 2. Materials and methods 2.1. The construction of T12-L1 model To establish a three-dimensional (3D) nonlinear complex FE model of T12-L1, 0.6-mm-thick computed tomography (CT) images of spine from a 50th percentile Chinese healthy female(no history of spinal disease) in the lying posture were obtained (Guo and Yin, 2019). The CT images were operated in the software (Mimics 10.0; Materialise Technologies, Leuven) to generate 3D geometry of TLJ. Then the geometry data was imported into CATIA V5R21 (Dassault Systèmes, VélizyVillacoublay) to create the 3D solid model of T12-L1. The model was meshed in pre-processing software of ANSA (BETA CAE Systems). The cortical bone and endplates modeled by 3-node shell elements consisted of the outer vertebral body and outer posterior elements. The interparts then were filled with cancellous bone core meshed by 4-node solid element. According to local variations of cortical and endplate thickness, the outer bone was subdivided into nine different thickness zones (Fazzalari et al., 2006; Roberts et al., 1997) (Fig. 1). Various locations
2.2. Material property All bony components were assumed to be homogenous and followed a strain rate-dependent elastoplastic Johnson-Cook material law (Table 1) that allowed computing von-Mises hardening with ductile damage until potential rupture (El-Rich et al., 2009). Before the equivalent stress reached to the yield stress, the material behaved as linear elastic. For higher stress values, the material behavior was plastic and the equivalent stress at constant temperature was calculated using the following Eq. (1):
σ = (A + Bε n)(1 + Clnε ∗̇ )
(1)
Fig. 1. The finite element model ofT12-L1 and Subdivision of the cortical bone, the cancellous bone into 9 components with different cortical thickness and into 7 components with different material properties. 30
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Table 1 Summary of the material properties of bone and IVD tissues used in the model (adapted from the reference (El-Rich et al., 2009; Lee et al., 2000; Schmidt et al., 2007; Wagnac et al., 2012)). Cancellous (A–G zone) and endplate centera
Material properties
Density (×10−6 kg/mm3) Modulus of elasticity, E (MPa) Poisson ratio, ν Yield stress, A (MPa) Harding modulus, B (MPa) Hardening exponent, n Failure plastic strain, εp Maximum stress (MPa) Strain rate coefficient, C Reference strain rate, ε0̇
Material properties
−6
A
B
C
D
E
F
G
1.8 93.7 0.25 1.95 8.5 1 0.082 2.65 0.533 0.008
1.8 93.7 0.25 1.95 7.0 1 0.06 2.3 0.533 0.008
2.0 93.7 0.25 1.95 8.5 1 0.082 2.65 0.533 0.008
2.0 93.7 0.25 1.95 8.1 1 0.08 2.6 0.533 0.008
2.5 93.7 0.25 1.95 12.5 1 0.104 3.25 0.533 0.008
2.5 93.7 0.25 1.95 12.5 1 0.104 3.25 0.533 0.008
1.8 93.7 0.25 1.95 7.0 1 0.06 2.3 0.533 0.008
Nucleus pulposus
3
Density (×10 kg/mm ) Poisson ratio, ν C10(MPa) C01(MPa) a
Cortical and endplate margin
0.2 4014 0.3 142.5 1338 1 0.071 190 0.272 0.008
Annulus groundsubstance
Low strain rate
High strain rate
Low strain rate
High strain rate
Fibers
1 0.495 0.12 0.03
1 0.495 31.8 8.0
1.2 0.45 0.18 0.045
1.2 0.45 11.8 2.9
Nonlinear
Bony endplate center material properties are the same as zone B.
Fig. 2. (A) Loading and boundary conditions of impact, (B) Loading and boundary conditions of quasi-static compressive test of the disc, (C) Compressive forcedisplacement curve of the disc, (D) loading and boundary conditions for validation of moment in a standard 3-D coordinate system, (E) degrees of left and right axial rotation, (F) degrees of flexion and extension, (G) Degrees of left and right lateral bending.
ε /̇ ε ̇ ε ̇ ≥ ε0̇ ε ∗̇ = ⎧ 0 ⎨ ⎩ 1 ε ̇ < ε0̇
Here, C10, C01 material constants; I1 and I2 the first and second invariants of the deviatoric component of the left Cauchy-Green strain tensor; J = V/V0, local volume ratio; d = 2/K, material incompressibility factor, K initial bulk modulus. But the law does not consider strain rate which has a significant influence on the mechanical behavior of IVD. Two sets of material parameters were adopted according to different strain rates (Schmidt et al., 2007; Wagnac et al., 2011). One set named as low strain rate was used in static or quasi static condition and the other set named high strain rate was used in high strain rate dynamic condition, such as the impact in this study. The property of collagen fibers was nonlinear based on a stressstrain curve (Shirazi-Adl et al., 1986). Since the stiffness and crosssection of fibers increase from the inner layers to the outer layers of the annulus, the revised ratio of the cross-sectional areas and elasticity
(2)
where, σ equivalent stress, A yield stress, B hardening modulus, n hardening exponent, ε ̇ current strain rate, ε0̇ reference strain rate and εpplastic strain (true strain). When εpreachesεmax, either for tension, compression or shear, shell elements would be deleted and the solid elements are not deleted but the deviatoric stress of solid elements are permanently set to zero. The IVD was governed by the hyper-elastic material law (Table 1) based on a first-order Mooney-Rivlin formulation with strain energy density equation, W is defined as follows.
W = C10 (I1 − 3) + C01 (I2 − 3) + (J − 1)/ d
(3)
K = 6(C10 + C01)/3(1 − 2ν)
(4) 31
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Fig. 3. The simulation results. (A)impact forces between the rigid ball and the rigid plate under different energy conditions, (B)the vertical displacements of centroid of the superior surface of T12 under different energy conditions, (C) Total contact forces between the facet surfaces under different energy conditions, (D) IDP changes under different energy conditions.
possible rebound (not impact tie contact). The ball was constrained and allowed only vertical motion to ensure vertical impact for minimizing the undesired flexion or extension. The inferior surface and spinous process of L1 were fixed in all degrees. The Radioss block100 explicit dynamic solver (Altair Hyperworks Inc. USA) was used to solve the investigation. The solution time was set to 6 ms under HE condition and 8 ms under LE and IME conditions and the running duration was about 4 h (Computer Configuration: Intel i52400, 3.10 Hz, 4GB RAM, Windows 7).
constants of the seven layers were adopted (Shirazi-Adl et al., 1986). The properties of ligaments were modeled in nonlinear material properties (Sharma et al., 1995). 2.3. Contact, loading and boundary conditions Due to the different mesh methods of bone, IVD and ligaments, a tied contact was created among them. Frictionless contact interfaces were used between adjacent facet joints to compute the contact forces. Frictionless contact interfaces were also created between different parts to avoid any possible penetration. All the loading and boundary conditions were set according to the literature (Kifune et al., 1995; Kifune et al., 1997; Panjabi et al., 1995; Panjabi et al., 1998; Qiu et al., 2006) (Fig. 2A). According to their experimental results, the TLJ initiated fracture in the impact energy of 2 kg ball (30 J) and BF in 4.3 kg (56 J). These were corresponded to the results that the energies replicating BF with no modification of the spinal dynamics are around 50 and 180 J (Germaneau et al., 2017; Jones et al., 2011; Panjabi et al., 2000). 1 kg, 2.3 kg and 4.3 kg were assigned to the rigid shell ball to represent low energy (LE), intermediate energy (IME) and high energy (HE) impact. In order to transfer the impact loads down to T12 uniformly and synchronously, a 3-mm thick rigid shell plate was created by the outline of T12 top surface. The shell plate was tied to the upper surface of T12. The rigid shell ball with the radius of 11 mm was placed over the rigid plate following its center vertically above the centroid of the plate with a distance of 5 mm between the shell plate and the lowest point of the ball. The ball was set at a 5.1 m/s downward velocity without gravitational acceleration to replace the experimental velocity of free falling from the height of 1.4 m for reducing computational cost. Contact was created between the rigid ball and rigid plate with
2.4. Validation of the model A quasi-static load was applied to the model of TLJ excluding posterior elements according to previous references (El-Rich et al., 2009; Zeng et al., 2013) (Fig. 2B). The loading rate was 1.26 mm/s and the loading duration was 1.2 s. KEREL order (kinetic energy relaxation) had been used to minimize the effect of kinetic energy and a curve of loaddisplacement was obtained. The curve was compared with the reported results (El-Rich et al., 2009; Zeng et al., 2013). The present model was validated under all physiological loading modes- flexion/extension, left/right lateral bending, and left/right axial rotation, by comparing these results with experimental results (Markolf, 1972; Oxland et al., 1992) and numerical results (Qiu et al., 2006) according to the range of motion (RoM) of the TLJ. The 7.5 Nm pure moment was applied to the TLJ in 3 incremental steps, 2.5 Nm, 5 Nm and 7.5 Nm, to create the relevant motion (Fig. 2D). 3. Results A compressive force-displacement curve of the 3-D nonlinear FE model of T12-L1was obtained by a velocity loading as showed in 32
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
initiated in vertebral foramen on the left-rear region on the superior cortical bone near the superior endplate of L1 at 1.85 ms. Then the fractures propagated along the circumferential direction. At 2.05 ms and 2.20 ms, fractures also occurred in the right-rear of cortical bone near the superior endplate and the posterior region of superior endplate of L1. The stress on cancellous bone adjacent to the fractural cortical bone was found to be zero meaning that the fracture occurred in the region. The fractural time was a little early than the cortical bone at 1.55 ms. Several fractural cancellous bone elements were located in the middle region of L1 on the left at the same time. The final fractural cancellous bone represented by dark blue elements was shown in Fig. 4C. In HE condition, more severe fracture occurred. The locations and the time of initial fractures were the same as those in IME, but the range of fracture was larger and the following propagation speed faster. At 2.8 ms, the fracture occurred in the middle of the antero-superior part of L1. With the rigid ball continuously compressing, the fractures in different locations connected together and formed a circle which divided L1 column into two parts. The stress distribution of cancellous bone was similar to that of IME, but the fracture was more severe. The superior and the anterior part of L1 collapsed toward the core of cancellous bone (Fig. 4D). As shown in Figs. 4 and 5, the cortical elements in the high stress region (red zone) exceeding the failure plastic strain would be deleted suddenly. The color of cancellous elements in the high stress region exceeding the failure plastic strain would turn to dark blue, but the elements keep deformation with zero stress output.
Fig. 2C. The curve fell within the quasi-static load-displacement experimental and FE corridors and showed a nonlinear compressive behavior. The predicted load-angular curves under various pure moments showed similar trends of incremental rotational motion and matched well with the published results in Fig. 2E–G. The validation study demonstrated that the present model was realistic. The impact forces between the rigid plate and the ball were gained (Fig. 3A). The HE reached to the maximum impact force of about 13KN. The minimum impact force was about 10.5 KN for LE. The vertical displacement was recorded by the centroid of the superior surface of T12 (Fig. 3B). According to Fig. 3A, the ball of LE first rebounded at about 1.8 ms and then the ball of IME rebounded at about 3.5 ms. The HE ball maintained a vertical motion without rebound. The maximum intradiscal pressure (IDP) caused by the three impacting energies were 44.5 MPa,50.2 MPa and 51.0 MPa respectively under the maximum impact forces (Fig. 3C). The rebound caused the IDP decreasing and vibrating in LE. The same case also happened in IME. The fractures in the cortical bone and endplate also decreased the IDP. The negative values of IDP represented the state of tension under impacting loads. The response from total contact force was markedly different under three impacting energies (Fig. 3D). In LE, the facet contact force exhibited a periodic motion, which was a cyclic contact between adjacent facets. The contact forces of IME and HE showed a vibrating attenuate tendency after the maximum impact force. The adjacent facets kept contact during the HE impacting time. The adjacent facets separated after the rebound and started a cyclic contact (Fig. 3D) under LE condition. In LE condition, the LE did not cause fracture of vertebral body. High stresses within cortical bone were located in the left and right rear region of L1 near the superior endplate margin and between the pedicles. The maximum stress reached to 237 MPa (Fig. 4A). High stress on cancellous bone was in the same region and the maximum stress was 11 MPa. In IME condition, the burst fracture occurred. Fractures were
4. Discussion This model used an accurate and practical method to investigate the mechanism of BF including the initiation, propagation and termination. This study not only defined high risk region owing to high-stress concentration (Lee et al., 2000; Qiu et al., 2006; Wilcox et al., 2004), but also reproduced the process of burst fracture. The load of vertical
Fig. 4. (A) The von Mises stress distribution in cortical bone under LE condition corresponding to the maximum vertical displacement of T12 (excluding the posterior part), (B) Initiation and propagation of cortical fracture in L1 under IME condition, (C) Initiation and termination of fracture of cancellous bone in L1 under IME condition, (D) Initiation and propagation of cortical fracture in L1 under HE condition. 33
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Fig. 5. The stress distribution of burst fracture of L1 cancellous bone under HE condition. (A) to (D) The fractures of cancellous bone at different time in different views, The time was set at 1.8 ms, 2.0 ms, 2.1 ms, 2.3 ms, 3.3 ms and 5 ms. (A) the left view, (B) the right view, (C) the front view, (D) the sagittal plane.
Depending on the failure strength of these regions, failure could be initiated in any of them. In this study, the initiation of fractures was on the posterior edge of superior endplate of L1 and adjacent cortical bone (including posterior cortical wall) (Fig. 3C and Fig. 4). The three results had shown varying degrees of damage. In the LE impacting condition, the TLJ was intact and no fracture was noted on the vertebral body. In the IME impacting (Fig. 4B–C), fractures occurred in both the cancellous bone and the cortical bone. The fractures initiated on the left posterior edge of superior endplate of L1 and the posterosuperior wall of vertebral body between the pedicle bases. Then the initial fractures trended to connect together and there was a new fracture on the right posterior edge of the superior endplate of L1. The fractures of cancellous bone adjacent to the fractured cortical bone occurred in the same time. A fractural zone in the middle of anterior column was observed where a potential fractured zone of burst fracture was described in the previous literature (Denis and Perry, 1982; Hakim and King, 1979; Panjabi et al., 1995). The fractures spread toward the core of the cancellous bone internally and would form an oblique line in the sagittal plane. Although the fractures occurred in IME, the vertebral body was a whole entity which had enough stiffness to rebound the impacting ball. The fragments of posterior wall invading the spinal cord recovered as before. There was no obvious fracture-dislocation of the TLJ. The structure of TLJ recovered after the impacting ball left the vertebral body. In HE (Fig. 5), the first half of fracture was similar to the fracture under the IME impacting condition. In the second half, with the impacting ball continuously compressing, the fracture continued to expand. Compared with the result of IME, a part of the anterior cortical bone and a part of the superior cancellous bone were isolated from the vertebral body along the fractural line in the sagittal plane resulting the
impacting ball was a transient variation with high to low velocity, which was closer to the fact in daily life but was different to the constant load (El-Rich et al., 2009; Fradet et al., 2014; Mustafy et al., 2014; Wagnac et al., 2011). The posterior ligament complex is seldom ruptured under axial compression unless the spine is subjected to forcible rotation or flexion-rotation (Holdsworth, 1970). The failure model of posterior ligament complex was simplified and was not analyzed in this study. Many researchers have offered classification and extensive description of fracture patterns in clinical observation. Magerl et al. (1994) described a total of 27 fracture patterns based on different single loads or combining forces. Vaccaro et al. (2005) evaluated the severity of spinal injury by measuring the parameters of spine. The FE method was a useful tool to verify the fracture patterns and provide the spinal injury severity score. Roaf (1960), Tran et al. (1995) and Qiu et al. (2006) suggested that increased pressure in the nucleus pulposus caused endplates to bulge toward the cancellous core resulting to crack which allowed the nucleus material to enter the vertebral body at a fast rate than the content could be expelled. The three results of this study appeared to contradict their hypotheses. The IDP only increased rapidly with the increased impact force. As the impact force decreased, the IDP in the nucleus never reached a maximum peak any longer. The bulge of IVD in the transverse direction was larger than it in the vertical direction meaning the rate of entering into the endplate was not great. This was consistent with the finding of Wilcox et al. (2004). The pressure and stress in the center of IVD were lower than those on the periphery which also supported the opinion. Lee et al. (2000) reported that high stresses could lead to fracture in the endplate and the posterior surface of the cortical wall.
Fig. 6. The patterns of burst fracture. (A)The FE result of this study, (B)Type B of burst fracture classified by Denis et al., (C)Type A of burst fracture classified by Magerl et al., (D)Lateral tomogram of a burst fracture. 34
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
and after burst fracture on the special model. The hypothesis that only anterior and middle columns fracture during BF in pure axial compression load was verified in this study. The burst fracture type under the HE condition was observed in clinical practice. The study suggested an integrated, detailed vision and insight into the mechanism, instability and treatment of TLJ burst fracture in clinical and experimental studies.
spinal cord stenosis. The vertebral body of L1 was divided into two parts (Fig. 6A). The same burst fracture was classified as type B by Denis and Perry (1982), which was described as disruption of superior endplate, loss of height, fracture of posterior wall of vertebral body (Fig. 6B and D)and was classified as type A3.1 by Magerl et al. (1994), which the upper or lower half of the vertebral body has burst while the other half remains intact (Fig. 6C). If the impact is not strong, the bone is within the elastic deformation, in ideal circumstances, that the bone will absorb the energy and start the reciprocating vibration, the ball can rebound. If the impact is strong, the plastic deformation will occur. The bone will absorb the energy and not release the energy which caused the plastic deformation. In this study, the ball rebound under LE and IME conditions. In LE condition, the velocity of rebound was 4.03 m/s, 5 J was absorbed by model. In IME condition, the velocity of rebound was 2.2 m/s, 9.8 J was absorbed by model. In HE condition, the plastic deformation was too large to converge in the analysis. Facet joints could play an important role during the impacting process. A hypothesis proposed by Wilcox et al. (2004) that the facet joints transmitted shear forces to the posterosuperior part of vertebral body. The fact that higher stress around the pedicle bases with high contact forces supported this hypothesis. The maximum IDP values in IME and HE were approximate which was different in LE. This may attribute to the sudden penetration of the nucleus into the gap of fractural bone, releasing the pressure. In this study, the trajectory of the ball is orthogonal to the plane of the endplate, so the rigid ball would not produce unexpected moment. The additional moment change the facet joint force and the morphology of IVD which have the contribution of different BF. This was discussed in Langrana et al. (2002) and had been completed in our following study. The numerical simulation has its limitations. It has been reported that the unique geometrical features characterizing each spine can determine significant changes in load distributions within the spine segments, so the lack of the diversity in resulting fracture patterns, various fractures under vertical impacting classified by the researchers (Denis and Perry, 1982; Holdsworth, 1970; Magerl et al., 1994) could not be obtained. This study could not compare between all the compressive fracture types. The study did not model the active and passive muscle tissues which may have influence on the TLJ kinematic responses. In high velocity dynamic loading, it was reasonable that there was no fluid flow in spinal segments as the aim of this study was focused on the mechanism of burst fracture under vertical impacting. It was different from the long-term (Argoubi and Shirazi-Adl, 1996; Schroeder et al., 2010; Heuer et al., 2007) compressive loads which needed visco-poroelastic property to study in a time frame of minutes or hours. But it still could be improved in a limited fluid-flow material property and independent mechanism under high velocities (Disilvestro and Suh, 2001; Oloyede and Broom, 1993). The load in this study was only axial compression and this could be improved by applying combined or complex loads to obtain more fracture patterns in clinical observations. The validation guaranteed the intended usage of FE model. The use of a verified model in static bending tests, especially, in quasi-static tests is one of the strengths of this study, but it still could be improved in dynamic validation which had a diverse response compered to in static condition. The L1-L5 dynamic experimental test results were obtained by Stokes and Gardner-Morse (2016), the TLJ segment experimental data was not found, and it should be completed in future work.
Declaration of Competing Interest The authors declare that they have no conflict of interest. Acknowledgement This work was supported by the National Natural Science Foundation of China (51875096, 51275082). References Alpantaki, K., Bano, A., Pasku, D., Mavrogenis, A.F., Papagelopoulos, P.J., Sapkas, G.S., Korres, D.S., Katonis, P., 2010. Thoracolumbar burst fractures: a systematic review of management. Orthopedics 33, 422–429. Argoubi, M., Shirazi-Adl, A., 1996. Poroelastic creep response analysis of a lumbar motion segment in compression. J. Biomech. 29, 1331–1339. Brandolini, N., Kapur, N., Hall, R.M., 2014. Dynamics of interpedicular widening in spinal burst fractures: an in vitro investigation. Spine J. 14, 2164–2171. Cotterill, P.C., Kostuik, J.P., Wilson, J.A., Fernie, G.R., Maki, B.E., 1987. Production of a reproducible spinal burst fracture for use in biomechanical testing. J. Orthop. Res. 5, 462–465. Daffner, R.H., Daffner, S.D., 2002. Vertebral injuries: detection and implications. Eur. J. Radiol. 42, 100–116. Dai, L.Y., Jiang, S.D., Wang, X.Y., Jiang, L.S., 2007. A review of the management of thoracolumbar burst fractures. Surg. Neurol. 67, 221–231. Denis, Francis, Perry, John F., 1982. The three-column spine concept and its significance in the classification of thoracolumbar spinal injuries. J. Trauma 22. Disilvestro, M.R., Suh, J.K., 2001. A cross-validation of the biphasic poroviscoelastic model of articular cartilage in unconfined compression, indentation, and confined compression. J. Biomech. 34, 519–525. El-Rich, M., Arnoux, P.J., Wagnac, E., Brunet, C., Aubin, C.E., 2009. Finite element investigation of the loading rate effect on the spinal load-sharing changes under impact conditions. J. Biomech. 42, 1252–1262. Fazzalari, N.L., Parkinson, I.H., Fogg, Q.A., Sutton-Smith, P., 2006. Antero-postero differences in cortical thickness and cortical porosity of T12 to L5 vertebral bodies. Joint Bone Spine 73, 293–297. Fradet, L., Petit, Y., Wagnac, E., Aubin, C.E., Arnoux, P.J., 2014. Biomechanics of thoracolumbar junction vertebral fractures from various kinematic conditions. Med. Biol. Eng. Comput. 52, 87–94. Fredrickson, B.E., Edwards, W.T., Rauschning, W., Bayley, J.C., Yuan, H.A., 1992. Vertebral burst fractures: an experimental, morphologic, and radiographic study. Spine 17, 1012–1021. Germaneau, A., Vendeuvre, T., Saget, M., Doumalin, P., Dupré, J.C., Brémand, F., Hesser, F., Brèque, C., Maxy, P., Roulaud, M., 2017. Development of an experimental model of burst fracture with damage characterization of the vertebral bodies under dynamic conditions. Clin. Biomech. 49, 139–144. Guo, L.X., Yin, J.Y., 2019. Finite element analysis and design of an interspinous device using topology optimization. Med. Biol. Eng. Comput. 57, 89–98. Hakim, N.S., King, A.I., 1979. A three dimensional finite element dynamic response analysis of a vertebra with experimental verification. J. Biomech. 12, 277–292. Heuer, F., Schmitt, H., Schmidt, H., Claes, L., Wilke, H.J., 2007. Creep associated changes in intervertebral disc bulging obtained with a laser scanning device. Clin. Biomech. 22, 737–744. Holdsworth, F.W., 1970. HoldsworthFFractures, dislocations, and fracture-dislocations of the spine. J Bone Joint Surg 52A: 1534-1551. J. Bone Joint Surg. 52, 1534–1551. Jones, H.L., Crawley, A.L., Noble, P.C., Schoenfeld, A.J., Weiner, B.K., 2011. A novel method for the reproducible production of thoracolumbar burst fractures in human cadaveric specimens. Spine J. 11, 447–451. Kifune, M., Panjabi, M.M., Arand, M., Liu, W., 1995. Fracture pattern and instability of thoracolumbar injuries. Eur. Spine J. 4, 98–103. Kifune, M., Panjabi, M.M., Liu, W., Arand, M., Vasavada, A., Oxland, T., 1997. Functional morphology of the spinal canal after endplate, wedge, and burst fractures. J. Spinal Disord. 10, 457–466. Langrana, N.A., Harten, R.D.R.D., Lin, D.C., Reiter, M.F., Lee, C.K., 2002. Acute thoracolumbar burst fractures: a new view of loading mechanisms. Spine 27, 498–508. Lee, C.K., Kim, Y.E., Lee, C.S., Hong, Y.M., Jung, J.M., Goel, V.K., 2000. Impact response of the intervertebral disc in a finite-element model. Spine 25, 2431–2439. Magerl, F., Aebi, M., Gertzbein, S.D., Harms, J., Nazarian, S., 1994. A comprehensive classification of thoracic and lumbar injuries. Eur. Spine J. 3, 184–201. Markolf, K.L., 1972. Deformation of the thoracolumbar intervertebral joints in response to external loads: a biomechanical study using autopsy material. J. Bone Joint Surg. (Am. Vol.) 54, 511–533.
5. Conclusion A detailed FE model of T12-L1 with nonlinear material property and failure criterion was built to research the mechanism of burst fracture under vertical impacting loads. The results provided the perspective on BF including the initiation, propagation and termination in detail, and the variations of stress, IDP, contact force and vertebral body before 35
Clinical Biomechanics 68 (2019) 29–36
L.-X. Guo and W.-J. Li
Pa 1976) 11, 914–927. Silva, M.J., Keaveny, T.M., Hayes, W.C., 1998. Computed tomography-based finite element analysis predicts failure loads and fracture patterns for vertebral sections. J. Orthop. Res. 16, 300–308. Stokes, I.A.F., Gardner-Morse, M., 2016. A database of lumbar spinal mechanical behavior for validation of spinal analytical models. J. Biomech. 49, 780–785. Trafton Jr., P.G., 1984. BC: computed tomography of thoracic and lumbar spine injuries. J. Trauma 24, 506–515. Tran, N.T., Watson, N.A., Tencer, A.F., Ching, R.P., Anderson, P.A., 1995. Mechanism of the burst fracture in the thoracolumbar spine. The effect of loading rate. Spine 20, 1984–1988. Vaccaro AR, Lehman RA, R John H, Anderson PA, Mitchel H, Rune H, James H, Marcel D, Kirkham W, Fehlings MG: A new classification of thoracolumbar injuries: the importance of injury morphology, the integrity of the posterior ligamentous complex, and neurologic status. Spine 2005, 30:2325. Wagnac, E., Arnoux, P.J., Garo, A., Elrich, M., Aubin, C.E., 2011. Calibration of hyperelastic material properties of the human lumbar intervertebral disc under fast dynamic compressive loads. J. Biomech. Eng. 133, 419–426. Wagnac, E., Arnoux, P.J., Garo, A., Aubin, C.E., 2012. Finite element analysis of the influence of loading rate on a model of the full lumbar spine under dynamic loading conditions. Med. Biol. Eng. Comput. 50, 903–915. White, A.A., Panjabi, M.M., 1990. Clinical Biomechanics of the Spine, 2nd edn. Lippincott, Philadelphia. Wilcox, R.K., Boerger, T.O., Allen, D.J., Barton, D.C., Limb, D., Dickson, R.A., Hall, R.M., 2003. A dynamic study of thoracolumbar burst fractures. J. Bone Joint Surg. Am. 85, 2184–2189. Wilcox, R.K., Allen, D.J., Hall, R.M., Limb, D., Barton, D.C., Dickson, R.A., 2004. A dynamic investigation of the burst fracture process using a combined experimental and finite element approach. Eur. Spine J. 13, 481–488. Willén, J., Lindahl, S., Irstam, L., Aldman, B., Nordwall, A., 1984. The thoracolumbar crush fracture. An experimental study on instant axial dynamic loading: the resulting fracture type and its stability. Spine 9, 624. Zeng, Z.L., Zhu, R., Li, S.Z., Yu, Y., Wang, J.J., Jia, Y.W., Chen, B., Cheng, L.M., 2013. Formative mechanism of intracanal fracture fragments in thoracolumbar burst fractures: a finite element study. Chin. Med. J. 126, 2852–2858. Zhao, F., Pollintine, P., 2009. Bd, Adams M, Dolan P: vertebral fractures usually affect the cranial endplate because it is thinner and supported by less-dense trabecular bone. Bone 44, 372–379.
Mcevoy, R.D., Bradford, D.S., 1985. The management of burst fractures of the thoracic and lumbar spine. Experience in 53 patients. Spine 10, 631–637. Mustafy, T., Elrich, M., Mesfar, W., Moglo, K., 2014. Investigation of impact loading rate effects on the ligamentous cervical spinal load-partitioning using finite element model of functional spinal unit C2-C3. J. Biomech. 47, 2891–2903. Natarajan, R.N., Andersson, G.B., 1999. The influence of lumbar disc height and crosssectional area on the mechanical response of the disc to physiologic loading. Spine 24, 1873–1881. Nicola, B., Nikil, K., Hall, R.M., 2014. Dynamics of interpedicular widening in spinal burst fractures: an in vitro investigation. Spine J. 14, 2164–2171. Oloyede, A., Broom, N.D., 1993. A physical model for the time-dependent deformation of articular cartilage. Connect. Tissue Res. 29, 251–261. Oxland, T.R., Lin, R.M., Panjabi, M.M., 1992. Three-dimensional mechanical properties of the thoracolumbar junction. J. Orthop. Res. 10, 573–580. Panjabi, M.M., Kifune, M., Wen, L., Arand, M., Oxland, T.R., Lin, R.M., Yoon, W.S., Vasavada, A., 1995. Dynamic canal encroachment during thoracolumbar burst fractures. J. Spinal Disord. 8, 39–48. Panjabi, M.M., Kifune, M., Liu, W., Arand, M., Vasavada, A., Oxland, T.R., 1998. Graded thoracolumbar spinal injuries: development of multidirectional instability. Eur. Spine J. 7, 332–339. Panjabi, M.M., Hoffman, H., Kato, Y., Cholewicki, J., 2000. Superiority of incremental trauma approach in experimental burst fracture studies. Clin. Biomech. 15, 73–78. Pintar, F.A., Yoganandan, N., Myers, T., Elhagediab Jr., A., 1992. SA: biomechanical properties of human lumbar spine ligaments. J. Biomech. 25, 1351–1356. Qiu, T.X., Tan, K.W., Lee, V.S., Teo, E.C., 2006. Investigation of thoracolumbar T12-L1 burst fracture mechanism using finite element method. Med. Eng. Phys. 28, 656–664. Roaf, R., 1960. A study of the mechanics of spinal injuries. Bone Joint J. 42-B, 810–823. Roberts, S., Mccall, I.W., Menage, J., Haddaway, M.J., Eisenstein, S.M., 1997. Does the thickness of the vertebral subchondral bone reflect the composition of the intervertebral disc? Eur. Spine J. 6, 385–389. Schmidt, H., Kettler, A., Heuer, F., Simon, U., Claes, L., Wilke, H.J., 2007. Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading. Spine 32, 748–755. Schroeder, Y., Huyghe, J.M., Donkelaar, C.C.V., Ito, K., 2010. A biochemical/biophysical 3D FE intervertebral disc model. Biomech. Model. Mechanobiol. 9, 641–650. Sharma, M., Langrana, N.A., Rodriguez, J., 1995. Role of ligaments and facets in lumbar spinal stability. Spine 20, 887–900. Shirazi-Adl, A., Ahmed, A.M., Shrivastava, S.C., 1986. Mechanical response of a lumbar motion segment in axial torque alone and combined with compression. Spine (Phila
36