Development and validation of lumbar spine finite element model

Introduction

The finite element method (FEM) has been widely used to improve road safety (John et al., 2022; Fang, Wang & Weggel, 2015; Sybilski & Małachowski, 2021; Bruski et al., 2019; Jones et al., 2016; Mendoza-Vazquez, Davidsson & Brolin, 2015). Simulations allow for investigating the behavior of the human body during a road accident. This feature offers the possibility of estimating the risk of fractures or damage to human soft tissues. One of the road traffic injuries with the most serious long-term consequences is a lumbar spine section injury (Zheng, Tang & Hu, 2018; Muller et al., 2014; El-Rich et al., 2009; Fradet et al., 2014; Eberlein, Holzapfel & Frohlich, 2004). Several new field studies have shown that lumbar spine fractures occur more frequently in new vehicle models than in previous models (Kaufman et al., 2013; Muller et al., 2014; Pintar et al., 2014; Doud et al., 2015). It is important to have reliable numerical models to explain this phenomenon.

Present crash test dummies were developed to assess injury criteria for the lower extremities, chest, neck, and head (Zheng, Tang & Hu, 2018). Their lumbar spine section arrangement was not based on actual human anatomy. Moreover, the validation for this section was not performed against any cadaver or human impact responses (Amiri, Naserkhaki & Parnianpour, 2020). Subsequently, despite the loading modes acting on the lumbar spine section being assessed in crash tests, they might not certainly follow the actual injury risk or fracture mechanism in frontal car crashes. Complete human body models, such as the Global Human Body Models Consortium (GHBMC), Virtual Vehicle Safety Assessment (ViVA), or Total Human Model for Safety (THUMS) are used to assess the risk of injury. Their advantage is a comprehensive approach to modeling, where the interaction of all body parts is taken into account. Moreover, they are usually validated against different experiments (Dreischarf et al., 2011; Guan et al., 2007; Heuer et al., 2008; Nachemson, 1960). The problem in their use may be the trade-off between the computation time and the level of complexity and discretization of a finite element (FE) model. This makes it problematic to study local phenomena by accurately identifying the location of failures. For this reason subsystems from a numerical model of the human body are used, such as models of the lumbar spine section. This allows for a more detailed analysis of the phenomena occurring in specific spine sections.

Various models of the lumbar spine can be found in the literature (Dreischarf et al., 2011; Xu et al., 2017; Schmidt et al., 2006; Renner et al., 2007; Wagnac et al., 2011). Before the model can be considered reliable, it is necessary to validate it against experimental results. Such validation has been described in previous works, such as in Zheng, Tang & Hu (2018); Östh et al. (2016); Wilke et al. (2001); Xu et al. (2017). Kiapour et al. (2012a) created a numerical model of the L3-S1 section based on computed tomography scans. This model has been validated using literature data for various load combinations (Goel et al., 2005; Kiapour & Goel, 2010; Kiapour et al., 2012b). The authors added a compression force to the model to simulate the weight of the upper body and muscles. A numerical model was used to analyze the movement of the lumbar spine after fusion surgery. Additionally, this model predicted the biomechanics of the spine after the insertion of an intervertebral disc implant. Zheng, Tang & Hu (2018) examined the lumbar spine injury mechanism in frontal car crashes. The effects of modification of thickness, stiffness, and angle of cushion, impact velocity, and coefficient of friction on lumbar spine injuries were analyzed. An article by Ayturk & Puttlitz (2011) presented several validation tests and indicated the need to extend the validation beyond the range of motion. Dreischarf et al. (2014) compared the results for several different numerical models that were considered validated. Thanks to the contribution of these works, it is possible to create a wider database of validated numerical models. This could potentially be used in the future to conduct statistical analyses, e.g., by examining the influence of various factors on the susceptibility of the spine to injuries. Due to the extensive biodiversity of geometry and material characteristics of the human body and still a limited amount of validation data, it is essential to fully share the results that can be obtained with the FEM.

In this study, a subsystem of the 50th percentile adult male (AM) THUMS v6.1 model (Toyota Motor Corporation & Toyota Central R and D Labs. Inc., 2021) was used, i.e., its lumbar spine section. Too stiff response of the THUMS lumbar spine model was observed in the flexion validation test for quasi-static and dynamic loading (see Fig. 4-44 and Fig. 4-46 in Toyota Motor Corporation & Toyota Central R and D Labs. Inc., 2021). Hence, the aim of the article is to improve the mechanical behaviour of the lumbar spine model in flexion and other validation tests. In order to achieve this goal nonlinear orthotropic material model was used for ligaments and modifications were introduced to the models of intervertebral discs to better reflect their anatomical structure and properties. Material parameters were taken from the literature and calibrated with respect to experimental flexion and compression tests, as they are of key importance in the analysis of road accidents (Richards et al., 2006; Ivancic, 2013; Pachocki et al., 2021; Tushak et al., 2023). In the present study, the set of validation tests were extended in relation the THUMS documentation (Toyota Motor Corporation & Toyota Central R and D Labs. Inc., 2021) by testing the range of motion in combined loading modes and testing the pressure level in the discs.

Materials and Methods

In this study, a three-dimensional (3D) FE model was used to simulate the biomechanics of the lumbar spine. The model was based on the 50th percentile AM THUMS v6.1, developed by Toyota Motor Corporation and Toyota R&D Labs Inc. (Toyota Motor Corporation & Toyota Central R and D Labs. Inc., 2021). The lumbar spine of the THUMS was extracted, modified, and subjected to various validation procedures. The modifications were applied both to the geometrical and material properties of the model. The FE mesh was refined; it consisted of 14,818 and 111,458 nodes that comprised 9,435 and 37,740 shell FEs, 49,972 and 503,712 solid FEs, and 11 and 17,478 beams for the THUMS and the current model, respectively. Nonlinear explicit dynamic analyses were performed for the entire system using LS-DYNA software (Hallquist, 2006; LSTC, 2017a; LSTC, 2017b). The general view of the model is depicted in Fig. 1A.

Boundary conditions and validation tests

Boundary conditions were applied to both models through load tables that are depicted in Fig. 1A. The bottom load table was fixed, whereas the specific force, moment, or prescribed motion was applied to the top table. To simulate epoxy fixing, some nodes of vertebra Th12 and L5 were constrained to the tables using a constrained nodal rigid body (*CONSTRAINED_NODAL_RIGID_BODY). The compression load was applied using cable beam elements (*SECTION_BEAM_TITLE, *MAT_CABLE_DISCRETE_BEAM), and the load path was assumed following the literature (Renner et al., 2007).

Table 2:

Thickness of ligaments of the lumbar spine model.

Ligament	Thickness L1–2, mm	Thickness L2–4, mm	Thickness L4–5, mm
Anterior Longitudinal Ligament (ALL)	2.5	3.5	1.5
Capsular Ligament (CL)	1.0	1.0	1.0
Ligamentum Flavum (LF)	5.0	5.0	5.0
Supraspinous Ligament (SSL)	4.0	3.0	2.0
Interspinous Ligament (ISL)	2.0	1.5	1.0
Intertransverse Ligament (ITL)	5.0	5.0	1.0
Posterior Longitudinal Ligament (PLL)	0.5	0.5	0.5

DOI: 10.7717/peerj.15805/table-2

Sixteen experimental tests were selected for comparison with numerical results. For 13 tests, the ranges of motion (ROMs) were compared between simulations and experiments; for the remaining three tests, the intradiscal pressures (IDPs) were compared. All experimental tests are summarized in Tables 3 and 4, where references to experimental data are given. The experimental results of IDPs were obtained for a single disc sample of L1-2 and L3-4 sections (Brinckmann & Grootenboer, 1991). For the IVDs from L2-3 and L4-5, the experiments were conducted on 6 samples. In the fourth column of Table 4 (Compression Force, N), the compressive forces adopted for the tests are listed. Table 4 summarizes the list of remaining high-speed tension and shearing tests from the work by Demetropoulos et al. (1998). The tests of numbers 1, and 2 were selected for the model calibration; the remaining tests were selected for the model validation. The model was considered valid when the range of motion of the whole lumbar spine section was within the range of the corresponding experimental data.

Table 3:

Quasi-static experimental tests for a range of motion (ROM) and interdiscal pressure (IDP).

Test number	Test name	Test type	Compression Force, N	Moment, Nm	Reference
1	Compression	ROM	1000	–	Renner et al. (2007)
10, 11, 12		IDP	300/1000/2000		Brinckmann & Grootenboer (1991)
2	Flexion	ROM	100	7.5	Xu et al. (2017)
3	Extension	ROM	100	7.5	Xu et al. (2017)
4	Lateral bending	ROM	100	7.5	Xu et al. (2017)
5	Torsion	ROM	100	7.5	Xu et al. (2017)
6	Flexion + Compression	ROM	1175	7.5	Xu et al. (2017)Renner et al. (2007)
7	Extension + Compression	ROM	500	7.5	Xu et al. (2017)Renner et al. (2007)
8	Lateral + Compression	ROM	700	7.8	Xu et al. (2017)Renner et al. (2007)
9	Torsion + Compression	ROM	720	5.5	Xu et al. (2017)Renner et al. (2007)

DOI: 10.7717/peerj.15805/table-3

Table 4:

Dynamic experimental tests for range of motion.

Test number	Test name	Test type	Maximum displacement, mm	Reference
13	Anterior shear	ROM	35	Demetropoulos et al. (1998)
14	Lateral shear		12
15	Posterior shear		35
16	Tension		2.5

DOI: 10.7717/peerj.15805/table-4

Results

This section presents the results of validation tests conducted for the two FE models: the current and THUMS v6.1. Both models were validated against experimental results obtained from the literature.

Range of motion in combined loading modes

The second group of tests consisted of four combined loading modes in which a compressive force was applied to the models simultaneously with a pure bending moment. In the first two tests, bending in the sagittal plane was analyzed. For a compression force of 1,175 N combined with a flexion moment of 7.5 Nm (see Fig. 5A), the rotations for the current model turned out to be consistent with the reference experimental values. Only the L3-4 segment failed 4.15° to fall within limits. In the case of the THUMS model, the rotations were outside the experimental data, where the differences ranged from 2.8° to 6.8°. For a compression force of 500 N combined with an extension moment of 7.5 Nm (see Fig. 5B), the current model fits very well with the experimental results. The THUMS model also achieved a good agreement with the results, only the response of L1-2 and L2-3 segments were too stiff, where the differences were 0.6° and 0.1°, respectively.

In the next performed tests, the combined loads acting in the coronal and transverse planes were analyzed. The results for these tests are summarized in Fig. 6. Both models successfully passed these tests, considering both the range of motion of the entire section and the ranges of motion for each functional spine unit (FSU).

Range of motion in dynamic tests

The last validation tests which were conducted are the examination of the ROM in the tension and shear tests. The responses of experimental data, the THUMS model and the current model for shearing tests are presented in Figs. 8A–8C. Both numerical models were close to or within the response corridor of shearing tests in the anterior and lateral directions. For the posterior shearing test, however, both of the models appeared to be too compliant compared to the experimental results. In all of the cases of shearing tests, the THUMS model was stiffer compared to the current FE model. In the 250 N tension test (see Fig. 8D), the responses of both models were inside the response corridor acquired from the experiment.

Discussion

In this article, a modified FE model of the lumbar spine based on the THUMS numerical model was validated against the experimental results of cadaveric tests available in the literature (Xu et al., 2017; Renner et al., 2007; Demetropoulos et al., 1998; Brinckmann & Grootenboer, 1991). The purpose of this study was to improve the response of the lumbar spine model in flexion by more realistic modelling of the soft tissues. An additional goal was to extend the validation of the THUMS model by testing the range of motion also in the combined loading modes and testing the pressure level in the discs. As a result, it is possible to investigate in detail deformation and the state of stresses in individual soft tissues. Numerical modeling of entire sections of the spine allows for analysis of risk injury that can occur in the human body during car accidents (Ivancic, 2013; Yoganandan et al., 2013; Pachocki et al., 2021; Kaufman et al., 2013).

Finite element models of the lumbar spine have been studied extensively using a variety of loading conditions, including compression loading. Compression loading can be used to study the behavior of the lumbar spine under loads that are similar to those typically encountered in daily activities. This can be important for understanding the effects of these loads on the spine and developing interventions to prevent or treat spinal disorders. A compressive force of 1,000 N was used for the calibration process, according to the experimental study. The current lumbar spine FE model presented a better response to this loading condition than the THUMS model.

Validation of lumbar spine FE models for tension and distraction loading conditions can be challenging. The lumbar spine is a complex structure, and these loading conditions can create complex stresses and strains within the spine that are difficult to measure experimentally. Additionally, the spine has a complex response under such loading, which is difficult to analyze in experimental tests fully. There have been several studies (Gay et al., 2008; Bae et al., 2020) that attempted to replicate the tension and distraction loading conditions on the lumbar spine using external loading devices, such as pull-tests on cadaveric spines, but these methods may also be insufficient to capture the true in-vivo response of the spine under tension/distraction. A literature review shows that there is a lack of tension validation for FE lumbar spine models (Xu et al., 2017). An experimental test by Demetropoulos et al. (1998) was done on cadaveric lumbar spine sections fixed at both ends and then tensioned. The same technique was performed for the current analysis. The simulation results for the modified model were consistent with the reported experimental data (within their range). Moreover, the obtained results were close to the average literature data. The THUMS model presented a stiffer response.

The predicted numerical results of the ROM analysis combining axial force and pure bending moments, matched with in vitro experimental data (Panjabi et al., 1994; Renner et al., 2007; Yoganandan et al., 1996). Additionally, a 1,175 N compressive force (Rohlmann et al., 2009) was used to simulate the model’s behavior during the carrying of external weight. Following in vivo experimental data from Rohlmann et al. (2001), a 100 N compressive follower force was added during the ROM in flexion, extension, and lateral bending motion. However, the compressive force leads to a significant ROM decrease for axial rotation. In extension, lateral bending, and axial rotation, the numerical results obtained for both numerical models were consistent with the experimental results. In all presented tests, the THUMS model predicted a slightly stiffer response. In addition, the lumbar spine response to complex load states also compares very well with the experimental ranges provided by the literature for the current model. However, the original model of THUMS v6.1 did not meet the condition of the range of experimental results for different loading modes of compression and flexion. In summary, these results motivated the authors to modify the THUMS model.

IDP validation was relatively more difficult due to the limited number of experimental studies available in the literature. Rohlmann et al. (2001) conducted in vitro research and determined the IDP in whole lumbar spine samples (L1–L5) by applying pure moments (axial rotation, extension, flexion, and lateral bending) plus a 280 N follower load. Nonetheless, the authors found that the IDP was smaller for the flexion compared to the extension for the same loading conditions, contradicting the results obtained by Wilke et al. (2008), Heuer et al. (2008) and Dreischarf et al. (2014) in the FE study of eight sections of the lumbar spine. In vitro experiments were also conducted on FSU samples, which contain one intervertebral disc and two vertebrae of the lumbar spine section (Heuer et al., 2008; Wilke et al., 2008). However, the IDP obtained for the whole section of the lumbar spine sample (Rohlmann et al., 2001) and for FSU specimens (Heuer et al., 2008; Wilke et al., 2008) differ significantly. For example, mean IDP in flexion (7.5 Nm) measured in L4-5 for the whole lumbar section (Rohlmann et al., 2001) and for FSU (Wilke et al., 2008) were 0.8 MPa and 0.4 MPa, respectively. The authors decided to validate the current FE model of the lumbar spine using results of in vitro experiments of discs conducted by Brinckmann & Grootenboer (1991). The current model predicted the intradiscal pressure in the NP very well for each functional spine unit and loading state. However, the pressure determined in the THUMS model was out of the range of the experimental data. This comparison indicates, that the current model improved the prediction of the experimental results. It is worth mentioning that the literature review indicates limited experimental data in the area of intradiscal pressure for NP.

In this study, the evaluation of the similarities metrics between experimental and numerical curves for both current and THUMS models were made. The study employed three statistical metrics. The results showed that in selected validation cases the current FE model performed better compared to the original THUMS v6.1 L-spine model in terms of selected similarity metrics, i.e., Pearson correlation, Sprague-Geers Comprehensive index, and WIFac rating. The result of a Pearson correlation indicates a better linear relationship with experimental data for the current model than the THUMS model. Sprague-Geers indices indicated that the results of the current model were closer to experimental data, in terms of magnitude and phase, compared to the THUMS model. Lastly, the WIFac showed that the overall shape of simulation curves in current model was closer to the expermintal results. By employing these three different metrics, the study provided a comprehensive analysis of the similarities between different sets of curves. Usage of these metrics allowed the evaluation of the relationships between curves with varying degrees of non-linearity and sampling frequency. Additionally, the study’s approach to averaging the metrics for each individual segment and subsequently for the whole section and corridors provides a more accurate representation of the similarities between the curves. Finally, the average value for all static and dynamic tests was shown. These results suggest that the current model may be more suitable for numerical analysis in biomedical research applications.

In numerical modeling, an important issue is the selection of a material law that best reflects the mechanical response of soft tissue. It is possible to describe the behavior of the structures in the lumbar spine section using various hyperelastic and other nonlinear material models such as Neo-Hookean, Mooney-Rivlin, or Yeoh models. Nevertheless, it is important to validate the kinematics of the spinal section as well as the mechanical response (pressure and forces) (Dreischarf et al., 2011). Dreischarf et al. (2014) presented different FE models subjected to pure compression and moments. They recommended analysis of the combined load of compression-torsion in experimental research. Hence, according to those recommendations, this article presents a wide range of loading modes acting on the lumbar spine, as well as an analysis of the pressure acting on the NP. In addition, a validation of tensile forces acting in the lumbar section was performed, which is not usually conducted. However, these loading modes are frequently observed during road accidents (Pachocki et al., 2021). The original contribution of the article is a separation of the NP from the surrounding bones and the AF. This has an impact on the response of this structure during different loading modes. The NP does not participate in the transmission of the tensile load, but in compression it does. Furthermore, the current model differs from the original THUMS model not only in constitutive laws and some values of material properties but also in the content and distribution of collagenous fibers in the AF ground substance and the anatomical thickness of ligaments. The described changes and modifications of the model resulted in a significant improvement in the model’s response to the given forces and moments during validation.

Numerical models are very useful in biomedical engineering, e.g., as a tool for the examination of risk analysis during car accidents. The differences in the results between the modified model and the original THUMS model show that one must be careful in drawing biomechanical conclusions from only one FE model for a particular population. The model needs to be validated to make more accurate analyses. Although the model may work well for examining the general response, it may not be suitable for detailed analyses. Nevertheless, the experimental analysis is still insufficient to explain the detailed failure mechanism, such as failure morphology and functions of soft tissues in the lumbar spine. Moreover, the anatomical and geometrical properties differ for experimentally tested samples (Mattucci et al., 2012), which may account for the differences in the obtained numerical results. The sensitivity of these parameters can influence the kinematic and mechanical response of the structure (Meijer et al., 2011; Robin, Skalli & Lavaste, 1994; Natarajan & Andersson, 1999; Lu, Hutton & Gharpuray, 1996; Niemeyer, Wilke & Schmidt, 2012). Moreover, the complex interactions and combinations of several geometrical and material parameters can govern the behavior of a model under particular combined loads. Hence, it is important to further develop comprehensive FE HBMs.

The presented work has its limitations:

• The 50th percentile of the AM THUMS v6.1 model was used in the analysis. Considering the biological biodiversity of the human population, the model yields results that properly describe the behavior of average male subjects. The current model can be used to investigate the biomechanics of a specific population. However, sensitivity analysis could provide importance for the development of numerical studies on the response of the lumbar spine section. It would show what structures, geometries, or material descriptions play a crucial role in the transmission of loads to the lumbar section.

• The presented validation was performed on the basis of the experimental data selected by the authors. However, more experimental data are available in the literature (Guan et al., 2007; Niosi et al., 2008; Pearcy, 1985; Yamamoto et al., 1989) that can be used in validation. The authors selected experimental studies conducted for a male subject due to the tested male model of the lumbar spine. However, due to the biodiversity of the population, it would be necessary to verify and validate the obtained data with other experimental works.

• A comprehensive analysis of the range of motion of the model in the sub-injurious range is an important first step in the validation of the model that may be used in crash analysis. However, it would be necessary to examine the model in conditions similar to those of a road accident, namely high strain rates, compression force, and moments (Minster, Lafon & Beillas, 2023; Yoganandan et al., 2023; Tushak et al., 2022; Tushak et al., 2023). The current validation set is not sufficient to fully evaluate the model’s ability to accurately predict the outcomes of a car crash.

Conclusions

A numerical model of the lumbar spine section was created based on the data available in the literature. The model passed most of the validation tests for various load combinations. The modifications introduced to the model resulted in a better agreement of the results with experimental data than for the original THUMS model. The separation of the NP from the AF and the surrounding bones changed the response of this structure to the action of compressive and tensile forces and allowed for obtaining the mechanical pressure in accordance with the literature. Moreover, the current model was validated against the kinematics of the spinal section, as well as the mechanical response. Considering the extensive complexity of this model, the level of this validation can be considered satisfactory. In the future, it is possible to expand the validation set with more extreme load cases.