A Numerical Investigation of Risk Factors Affecting Lumbar Spine Injuries Using a Detailed Lumbar Model

Recent field data showed that lumbar spine fractures occurred more frequently in late model vehicles than early ones in frontal crashes. However, the lumbar spine designs of the current crash test dummies are not accurate in human anatomy and have not been validated against any human/cadaver impact responses. In addition, the lumbar spines of finite element (FE) human models, including GHBMC and THUMS, have never been validated previously against cadaver tests. Therefore, this study developed a detailed FE lumbar spine model and validated it against cadaveric tests. To investigate the mechanism of lumbar spine injury in frontal crashes, effects of changing the coefficient of friction (COF), impact velocity, cushion thickness and stiffness, and cushion angle on the risk of lumbar spine injuries were analyzed based on a Taguchi array of design of experiments. The results showed that impact velocity is the most important factor in determining the risk of lumbar spine fracture (P = 0.009). After controlling the impact velocity, increases in the cushion thickness can effectively reduce the risk of lumbar spine fracture (P = 0.039).


Introduction
Safety designs in newer cars generally provided better protections to occupants than those in older cars. However, several recent field data analyses have shown that lumbar spine fractures occurred more frequently in late model vehicles than the early ones in frontal crashes [1,2].
This increasing trend of lumbar spine fractures in frontal crashes is extremely concerning, because none of the current crash test programs, including FMVSS 208, US-NCAP, and IIHS, considered lumbar spine injury in their safety evaluation process. This is partially due to the fact that none of the current crash test dummies can accurately estimate lumbar spine fracture risks. Current crash test dummies were designed to focus on estimating injury risks to the head, neck, chest, and lower extremities, and their lumbar spine designs were not based on real human anatomy and were not validated against any human/cadaver impact responses. As a result, although lumbar spine loadings could be measured in a crash test, they may not necessarily reflect the real injury mechanism or injury risk in frontal crashes. In addition, lumbar spines of the available FE human models, such as GHBMC (Global Human Body Model Consortium) model and THUMS (Total Human Model for Safety) model, have never been validated previously against cadaver tests. More recently, Arun et al. [3] applied the stiffness values obtained from cadaver tests directly to the lumbar spine of the GHBMC simplified model. However, it is a rigid bodybased lumbar spine model. It is lacking of detailed anatomical structures of the lumbar spine and cannot estimate strain and stress in the vertebrae.
Factors affecting the lumbar spine fracture have been discussed extensively in the literature. Results indicated that older occupants (65+ years) were five times more likely to sustain spinal injury compared to younger occupants [4], women had more lumbar spine fractures than men [5], and lower bone quality is associated with an increased number of lumbar spine fractures within the CIREN cases analyzed [6]. Lumbar spine posture was also found to be an important factor affecting lumbar spine injuries. It was reported that slouched posture may cause an increase in stress on the lumbar vertebrae [7] and more reclined postures are associated with a higher lumbar vertebrae fracture risk by reconstructing real-world motor vehicle crashes [8]. Crash pulse magnitude and shape are also crucial for determining lumbar spine injury risk. For example, cadaver tests under axial loading showed that triangular pulses produced pelvic fractures with no lumbar spine fractures, while a sigmoid-shaped pulse produced no pelvic fractures but did produce an L1 burst fracture [9]. In addition, more severe crash pulses may lead to higher lumbar spine injury risks [10].
The objectives of the present study were to develop a detailed FE lumbar spine model, validate it against cadaveric tests, and use it to investigate the effects of changing the coefficient of friction (COF), impact velocity, cushion thickness and stiffness, and cushion angle on the risk of lumbar spine injuries. This study could provide better understanding of how to design countermeasures to reduce occupant lumbar spine injuries in the new generation of vehicle models.

Lumbar Spine Model Development and Validation
2.1. Model Development. In this study, the GHBMC model was selected as the baseline model to be modified due to its accurate geometrical representation. In the original GHBMC lumbar spine model (Figure 1(a)), vertebrae (T12-L5) were represented by a rigid material, intervertebral discs were modeled by shell elements (Figure 2(a)), and no ligaments were modeled. In the modified model ( Figure 1(b)), the intervertebral disc was separated into three parts including nucleus, annulus, and collagenous fibers (Figure 2(b)). A 0.5 mm layer of shell element was added to the outside of the original vertebrae, representing the cortical bone. Detailed ligaments modeled by beam elements were added between the vertebrae. Nucleus, annulus, and cancellous bones were modeled by solid elements. Collagenous fibers and ligaments were modeled by beam elements. The facet joint articulations were simulated as the tiebreak contact between adjacent endplates. Note that muscles were ignored in the current study. The modified lumbar model (Figure 1(b)) was composed of 64,338 elements and 12,295 nodes. The cancellous bones and cortical bones were modeled by MAT_SIMPLIFIED_-JOHNSON_COOK in LS_DYNA [11,12] with 1.3% and 0.94% effective plastic failure strains, respectively. The endplates of the discs were modeled by an elastic material, and   the endplates of vertebrae were modeled by MAT_SIMPLI-FIED_JOHNSON_COOK. The nucleus and annulus were modeled by MAT_MOONEY-RIVLIN_RUBBER [11,12]. Accordingly, the shear modulus constant A of the nucleus and annulus were 0.64 MPa and 0.24 MPa, and the constant B of nucleus and annulus were −0.16 MPa and −0.06 MPa, respectively [12]. A linear elastic material was used for ligaments. Selections of stiffness of ligaments (Table 1) were based on a previous study [13]. The collagenous fibers were represented by a nonlinear load-displacement curve obtained from the literature [14]. Because external lamellae are stiffer than internal lamellae, the fibers in different layers were weighted according to a previous research [15]. To tune the model responses, the material properties of the bones and ligaments were adjusted slightly to match the test results during model validations. The material properties of all the parts in the modified lumbar model are listed in Table 2.

Validation against Nonfailure Tests. Cadaveric tests from
Demetropoulos et al. [17] were used to validate the modified lumbar model. It is important to adjust the initial posture of the lumbar spine to be consistent with the experimental data before validation [18]. Therefore, a presimulation was conducted to adjust the spine curvature to match the specimen pretest condition, as shown in Figure 3. Impact simulations were performed at different loading modes corresponding to the testing conditions. As a result, six simulations under compression, anterior shear, posterior shear, lateral shear, extension, and lateral bending were performed with the modified lumbar spine model. Each simulation was run at a displacement rate of 100 mm/sec. The maximum displacements were set to be the same as those in the tests. During these simulations, T12 was fixed to the fixture on the top, and L5 was attached to the fixture at the bottom.
In the compression and shear conditions, load was applied from the lower fixture with the upper fixture rigidly constrained. In the bending condition, displacement was applied from the lower fixture, and a bending moment was applied to the superior end of the lumbar spine by a cable. Force, moment, and angular displacement were measured at the upper fixture. Figure 4 shows load-displacement curves in each loading condition for the modified lumbar spine model. In compression, posterior shear, lateral shear, extension and lateral bending, simulation results shown in red line fell within test corridors reasonably well.

Validation against Failure
Tests. Cadaver tests with tissue failure conducted by Duma et al. [19] were used to validate the modified lumbar model. The validation simulations were set up under the same configurations as those in the dynamic compression tests ( Figure 5). The first simulation used a lumbar spine motion segment ( Figure 5(a)), and the second simulation used the entire lumbar spine ( Figure 5(b)). During the simulations, the superior end of the lumbar spine was attached to the upper plate, and the inferior end was attached to the lower plate. A prescribed motion was applied on the upper plate while the lower plate was fixed. All the failure simulations were loaded at 1.0 m/s. Force and moment were calculated in these simulations under different loadings.
As shown in Figure 6, the model-predicted forcedisplacement curve for the motion segment configuration fell within the experimental corridor. Good correlations between the simulation and the test for the entire lumbar spine   compression were also achieved, as shown in Figure 7. In addition, the locations of fractures in the simulation were consistent with those in the test (underneath the T12 and L3), as shown in Figures 7(c) and 7(d).

Design of Experiment Analysis
A DoE analysis based on Taguchi Array was used to study the effects of multiple factors on the risk of lumbar spine injuries.   In this study, 5 factors were investigated, including cushion stiffness, cushion thickness, cushion angle, the coefficient of friction, and impact velocity. The setup of the simulations in the DoE analysis is shown in Figure 8 [20]. The upper mass was attached to the upper platform to simulate the mass of the torso, head-neck, and upper extremities. The interaction between the upper platform and the upper fixture was in a form of a laterally oriented cylinder. The lumbar spine was fixed at cranial and caudal ends to the upper fixture and lower fixture, respectively [20]. A cushion foam, of which the thickness was defined as h, was attached to the fixed plate as shown in Figure 8, and its density was varied for the DoE analysis. The fixed plate was constrained at a cushion angle of θ. The initial impact velocity was applied to the upper mass, upper platform, and impact cylinder, through the upper fixture, lumbar spine specimen, and lower fixture, and finally reached the lower platform, cushion foam, and the fixed plate ( Figure 8).
The initial setup conditions of the impact velocity, cushion angle, thickness, and density were 0.1 m/s, 0 degree, 20 mm, and 20 kg/m 3 , respectively. The coefficient of friction (COF) between the lower platform and the seat cushion was ranged from 0 to 0.9. Four different levels were assigned for each factor (Table 3). Because most lumbar injuries were reported as vertebral fractures, the maximum principal strain in the bony structures was selected to evaluate the risk of lumbar fracture [21]. These setups resulted in a total of 16 simulations based on the Taguchi Array, as shown in Table 4. One-way analysis of variance (ANOVA) and analysis of covariance (ANCONA) were performed using SPSS 20.0.
The average values of the maximum principal strains of different factor categories are shown in Figure 9. Impact velocity (P = 0 009) had the most significant influence on the risk of lumbar injuries. With impact velocity varying between 0.1 and 0.5 m/s, the maximum principal strain fluctuated between 0.0032 and 0.0038, but it increased rapidly to 0.0053 when impact velocity rose to 0.7 m/s. Except the impact velocity, no other factor was statistically significant (P > 0 05).
Because the impact velocity dominated the lumbar fracture risk in the Taguchi array, a one-way ANOVA was then conducted for the remaining 4 factors by controlling the impact velocity as a constant. The average values of the maximum principal strains of different factor categories are shown in Figure 10. The cushion thickness (P = 0 039) became significant for the lumbar injuries. An increase in the cushion thickness will lead to more energy absorption and in turn lower lumbar spine fracture risk. Even though other factors were not significant, some general trends have also been noted. For example, with an increase in the COF, the lumbar spine fracture risk generally increased. For the cushion angle, the highest injury risk occurred when the lumbar spine orientation was perpendicular to the cushion. These trends are widely consistent to the findings from other studies on lumbar spine injuries [22] and cervical spine injuries [23].

Conclusions
A detailed modified FE lumbar spine model based on the GHBMC model was developed and validated against available cadaveric tests which appeared in the literature. Risk factors affecting lumbar spine injuries were investigated using the modified model. Results of the DoE analysis demonstrated that the impact velocity is a significant factor influencing the lumbar injuries. After controlling the impact velocity, cushion thickness is another significant factor influencing lumbar injuries. An increase of cushion thickness or decrease of cushion stiffness will reduce the lumbar spine fracture risk. Additionally, minimizing the COF between the padding and the lumbar spine can reduce the lumbar spine injury risk.
One of the limitations of this study is that human factors such as BMI (body mass index), sex, and age were not considered. Older occupants have higher risk of lumbar spine fractures, and their lumbar spine curvatures may be different to the younger adults [24]. In addition, muscles were not included in the current lumbar model. Future studies may investigate the effects of human factors and muscle activations on the lumbar spine injury risks using parametric human models [25,26].

Conflicts of Interest
The authors declare that they have no conflicts of interest.