Modular use of human body models of varying levels of complexity: Validation of head kinematics

ABSTRACT

Objective: The significant computational resources required to execute detailed human body finite-element models has motivated the development of faster running, simplified models (e.g., GHBMC M50-OS). Previous studies have demonstrated the ability to modularly incorporate the validated GHBMC M50-O brain model into the simplified model (GHBMC M50-OS+B), which allows for localized analysis of the brain in a fraction of the computation time required for the detailed model. The objective of this study is to validate the head and neck kinematics of the GHBMC M50-O and M50-OS (detailed and simplified versions of the same model) against human volunteer test data in frontal and lateral loading. Furthermore, the effect of modular insertion of the detailed brain model into the M50-OS is quantified.

Methods: Data from the Navy Biodynamics Laboratory (NBDL) human volunteer studies, including a 15g frontal, 8g frontal, and 7g lateral impact, were reconstructed and simulated using LS-DYNA. A five-point restraint system was used for all simulations, and initial positions of the models were matched with volunteer data using settling and positioning techniques. Both the frontal and lateral simulations were run with the M50-O, M50-OS, and M50-OS+B with active musculature for a total of nine runs.

Results: Normalized run times for the various models used in this study were 8.4 min/ms for the M50-O, 0.26 min/ms for the M50-OS, and 0.97 min/ms for the M50-OS+B, a 32- and 9-fold reduction in run time, respectively. Corridors were reanalyzed for head and T1 kinematics from the NBDL studies. Qualitative evaluation of head rotational accelerations and linear resultant acceleration, as well as linear resultant T1 acceleration, showed reasonable results between all models and the experimental data. Objective evaluation of the results for head center of gravity (CG) accelerations was completed via ISO TS 18571, and indicated scores of 0.673 (M50-O), 0.638 (M50-OS), and 0.656 (M50-OS+B) for the 15g frontal impact. Scores at lower g levels yielded similar results, 0.667 (M50-O), 0.675 (M50-OS), and 0.710 (M50-OS+B) for the 8g frontal impact. The 7g lateral simulations also compared fairly with an average ISO score of 0.565 for the M50-O, 0.634 for the M50-OS, and 0.606 for the M50-OS+B. The three HBMs experienced similar head and neck motion in the frontal simulations, but the M50-O predicted significantly greater head rotation in the lateral simulation.

Conclusion: The greatest departure from the detailed occupant models were noted in lateral flexion, potentially indicating the need for further study. Precise modeling of the belt system however was limited by available data. A sensitivity study of these parameters in the frontal condition showed that belt slack and muscle activation have a modest effect on the ISO score. The reduction in computation time of the M50-OS+B reduces the burden of high computational requirements when handling detailed HBMs. Future work will focus on harmonizing the lateral head response of the models and studying localized injury criteria within the brain from the M50-O and M50-OS+B.

KEYWORDS:

• Biomechanics
• computational modeling
• human modeling
• injury

Introduction

Computational human body models (HBMs) can be used as tools to study injury prediction in accelerative loading environments ranging from vehicular crash (Arun et al. 2016), to aerospace applications (White et al. 2014), to sports applications (Post et al. 2014). The Global Human Body Model Consortium (GHBMC) is one such model currently utilized for vehicle crash response. The GHBMC is an industry-sponsored and government-supported consortium with the aim of developing computational human models for the blunt injury environment.

The GHBMC average male occupant (M50-O v4.4) finite-element model was used in this study. The development of this HBM was based on a multimodality medical image and external anthropometry data set of a volunteer representing a 50th percentile male in terms of height (174.9 cm) and weight (78.6 ± 0.77 kg), described by Gayzik et al. (Gayzik et al. 2011; Gayzik et al. 2012). The model has since undergone numerous validation simulations at both the regional (e.g., DeWit et al. 2012; Li et al. 2010; Shin et al. 2012; Soni et al. 2015) and full body levels (e.g., Hayes et al. 2014; Toyota 2010; Yang et al. 2006). Additional information on the development of the model can be found in the GHBMC M50 user's manual (GHBMC 2011).

The complexity of the detailed M50-O model requires significant computational resources, which motivated the development of a simplified average male occupant model (GHBMC M50-OS). The M50-OS yields similar response characteristics as the detailed model with a drastically reduced computational cost. The M50-OS retains the same habitus and rigid bone structures as the detailed model but has reduced mesh density and greatly simplified soft-tissue structures. For example, the organs are homogenized and the musculature is either represented as one-dimensional beam elements or considered part of the soft-tissue envelope (Schwartz et al. 2015). The neck model of the M50-OS is based on work by Dibb (2011). Research has shown that the M50-OS model runs roughly 32 times faster than the M50-O (Decker et al. 2016).

While details of the simplified model can be found in Schwartz et al. (2015), the present study focuses on head kinematics; thus, the brain approach is revisited briefly. Rather than the detailed M50-O brain (Mao et al. 2013), the M50-OS has a mass node of 1.6 kg at the center of gravity (CG) of the brain, constrained to the skull. Previous studies have demonstrated the ability to modularly incorporate the validated GHBMC M50-O brain model into the simplified model (GHBMC M50-OS+B). Compute times were reduced by a factor of 9 compared to the time required to for the detailed model (Decker et al. 2016). Despite demonstrating this capability, these models have yet to be independently validated.

The Naval Biodynamics Laboratory (NBDL) conducted a study on head–neck response of human volunteers in frontal and lateral loading conditions. This test program consisted of 16 subjects with a test database of 72 lateral and 119 frontal sled tests (Wismans et al. 1986). The NBDL data set is a widely used kinematic data set for model validation. While validation using these studies have been conducted on an isolated head and neck structure as part of the neck model validation of the GHBMC M50-O model, the accelerations were applied directly at the first thoracic vertebral body (T1) (GHBMC 2011). The objective of this study therefore is to validate the head and neck kinematics of the GHBMC M50-O and M50-OS (detailed and simplified versions of the same model) against NBDL human volunteer test data in frontal and lateral loading using full-body simulations. Given unknowns about the precise belt system for these tests, a sensitivity study on a number of parameters is provided for the frontal configuration. Furthermore, the effect of modular insertion of the detailed brain model into the M50-OS is quantified.

Methods

The overall approach for the study was to compare data from published human subject tests at 15g frontal, 8g frontal, and 7g lateral against the M50-O, M50-OS, and M50-OS+B simulations. The methods section is broken into two parts. First, we present data procurement and preparation methods, followed by the computational model preparation and simulation methods.

Data procurement and preparation

An anthropometric comparison of human volunteer mean from the NBDL tests and each HBM was conducted (Table A2, see online supplement). Volunteer subject head masses were determined using published regression equations (McConville et al. 1980) and were reported by Wismans et al. (1986). Head center of gravity of the volunteers in the study was defined by Beier et al. and the anatomical origin of T1 was defined as the anterior superior corner of the first thoracic vertebral body (Beier et al. 1979; Thunnissen 1995). Model data were measured from the models themselves, but head mass was measured following the same section plane approach as McConville (Vavalle et al. 2013). The reported location of the head CG is the location of the head CG node in each of the models, which was determined computationally.

The three-dimensional kinematics of the head and first thoracic vertebral body (T1) were experimentally monitored in the study by anatomically mounted clusters of accelerometers and photographic targets (Thunnissen 1995; Wismans et al. 1986). Mean sled accelerations and volunteer kinematics were obtained for each of the three test conditions from the online National Highway Traffic Safety Administration (NHTSA) Biomechanics Test Database (NHTSA 2012). Representative (average) curves and corridors presented are the results of a reanalysis of all data (sled acceleration, human body accelerations) using methods described by Gayzik et al. (2015), which is work based on Nusholtz et al. (2013).

Model preparation

NBDL test positioning and belting

All models were run on the Distributed Environment for Academic Computing (DEAC) high-performance computational cluster at Wake Forest University using 48 central processing units (CPUs) and LS-DYNA R.7.1.2 (LSTC 2014; Livermore, CA). The setup for each respective NBDL test was reconstructed to match occupant positioning and the belting setup with available data provided in the experimental description. In the frontal tests, the human volunteers were seated in an upright position with a leg angle of 134 to 140 degrees. The occupants were restrained by a five-point belting system, with two shoulder straps connecting to a lap belt and crotch strap attached to the underside of the back rest. A 3-inch-wide Dacron belt was used for the belting system, with wrist and foot restraints used in all testing. The same setup and belting system was used for the lateral tests, with an additional 25-cm-wide shoulder strap and a lightly padded wooden board placed against the right shoulder to limit lateral motion of the subject. The test buck and belts were modeled as described in the following.

Significant effort was made to match the initial position of the models to the volunteer position (to the extent possible based on experimental description), as this has a high influence in resulting kinematics. Exact positional matching for the whole body was not possible due to a lack of absolute measurements of bony landmarks on the NBDL volunteers. However, the spatial relationship between the head CG and T1 was provided by the NBDL study. It was deemed appropriate to match HBM position to the volunteers based on seating angle, leg angle, and head-to-T1 relationship.

Human model positioning

Positional matching was accomplished for each HBM with a settling simulation that addressed each of these target areas. The seated position was obtained by rotating the seatback angle of a generic setup to 90 degrees, while the legs are simultaneously pushed anterior–posterior by a knee bolster and footrest 70 mm in order to keep the back flush against the seatback and retain the proper leg angle. Gravity was active during the settling simulation. Proper positional matching of the head and T1 was accomplished using the keyword Boundary_Prescribed_Final_Geometry, where nodes can be displaced along a straight line trajectory from their initial positions to the final location defined in the input (LSTC 2014). The resulting head-to-T1 relationship after positioning was compared between the three models and to the volunteer data (Table A1, see online supplement). The angle between the head CG and T1 was matched between models and volunteer data, but lengths were not explicitly matched as no scaling of the models was done. The neck length of the model and reported lengths of the subjects are in Table A1 (see online supplement). The 7g lateral setup required additional positioning to reposition the M50 models' arms inward and to match the experimental setup by limiting lateral motion of the body during the simulation pulse. Because no head positioning data for lateral were available we assumed that position would be identical to the frontal setups. The setup includes a lateral plate pressed against the right shoulder of the HBM similar to that in the NBDL testing. The plate extends down past the hip of the HBM to prevent lateral pelvis motion and improve model stability. The exact dimensions of the lateral board were not released by the study. An additional plate on the left shoulder is used to symmetrically compress the arms to the HBMs and has no effect on the HBM response during the simulation. Arm compression was accomplished with an inward displacement of 45 mm by both plates during the settling simulation.

Human model belting

The models were belted using a mixed one-dimensional (1D) and two-dimensional (2D) element setup. The width and thickness of the 2D belt elements were matched to the belt used in the NBDL study. The material property for the 2D belt was set as MAT_FABRIC and the 1D belt as MAT_SEATBELT. The belts used in the model use empirical data obtained from a validated belted simulation (Davis et al. 2016) and yield ∼9 kN of force at 6% strain. These seatbelt response values are consistent with data provided by original equipment manufacturer (OEM) collaborators through the GHBMC. The 1D seatbelt elements were implemented at the edges of the belt to allow for implementation of belt pretensioner and retractor functions. These were not meant to match the experimental setup per se, but rather were used to tighten the HBMs against the seat. A slip-ring system was developed for the 1D seatbelt elements of the crotch strap to guide the belt under the occupant during pretensioning.

Boundary conditions

Following belting, the models were run in the 8g frontal, 15g frontal, and 7g lateral simulations. All parts of the buck were defined as rigid parts. The corresponding sled velocity pulse determined in the data preparation phase was applied to the buck using the Boundary_Prescribed_Motion card. The frontal and lateral simulations were run with the M50-O, M50-OS, and M50-OS+B with neck muscle activation to approximate the living volunteer response. The same curve used for active musculature in the GHBMC M50-O was used in the M50-OS and M50-OS+B with the curve beginning at 134 ms of the simulation (GHBMC 2015; Siegmund et al. 2003). This includes a 74-ms response time observed in Siegmund et al.

Several details were not specified related to the belt system, such as belt loads, seat reaction loads, belt pretension values, and belt webbing force versus displacement data. To model the belt system as realistically as possible, some additional settings were used in the computational model. A modest amount of tension in the belts at the time of the pulse was used to represent a taught belt. While we understand the belt system used in the actual experiments did not employ automatic pretensioners, these were used to attempt to couple the body as closely as possible to the buck during the simulation. Retractors were implemented as part of the LS-DYNA seatbelt approach and to account for some slack that was likely inherent in the system. Reasonable values of peak pretension were used in the simulations for the shoulder (150 N), lap (400 N), and crotch (400 N). Peak values were chosen to be high enough to hold the HBM in place while staying within realistic values. A sensitivity study was conducted to attempt to quantify the effect of these assumptions. Variables observed to have the highest impact on head kinematics in preliminary simulations were chosen for analysis. The parameters included were pretensioning time, belt stiffness, the time at which any slack in the system was locked, and time of neck muscle activation.

The time over which these systems must operate were added through reasonable assumptions. As the time of maximum pretension, 65 ms was selected for the baseline simulations. Velocity of the buck at this point was observed to be approximately 10% of the maximum buck velocity and was considered to be an acceptable time for maximum pretension. The effect of maximum pretension time was tested from 50 ms (∼2% maximum buck velocity) to 65 ms in increments of 5 ms. The effect of the loading and unloading characteristics of the belt was tested with a sensitivity analysis of ±10% and 20% scaling factors on the loading and unloading ordinate from the experiments, which are directly implemented in the model. Sensitivity analysis of neck muscle activation involved offsetting the current muscle activation curve to represent three conditions: the default condition described already, an early onset of bracing, and no bracing. Early bracing removed the 74-ms response time, and no bracing was modeled without muscle activation. The inherent slack in the belt system was tested by varying the shoulder-belt retractor locking time. This varied from 65 ms (no slack) to 120 ms, with 100 ms as the baseline. Visual comparison of the sensitivity analysis can be seen in Figure A7 (see online supplement). This analysis was conducted for a single case in which the belt system was thought to have the most effect, 15g frontal.

The targeted kinematic measures for this study include angular head CG acceleration, linear resultant head CG acceleration, and linear resultant T1 acceleration. Angular head acceleration is defined as +Y for frontal and +X for lateral motion per the SAE J211 standard. Head accelerations were taken from a node at the head CG for which the location and local coordinate system are interpolated based on head motion. T1 acceleration data were obtained from the rigid body output (RBDout). Each of these respective measures was output from each model during the 8g frontal, 15g frontal, and 7g lateral simulation. For each test setup, the model results were compared to target corridors developed from reanalysis of the NBDL data. Anthropometric locations used to measure kinematics were cross-checked between the HBMs and the definitions given by the NBDL study.

Quantitative comparison

The quantitative comparison of results was conducted using the ISO TS 18571 approach (Barbat et al. 2013). There are three submetrics (phase, magnitude, and slope) and a corridor score used to calculate ISO, all ranging from 0 to 1 (perfect score). The submetrics represent how well a curve matches different aspects regarding the shape of the volunteer mean curve. The corridor score is calculated using the algorithm from the CORA objective evaluation technique and compares the model curve to a set of inner and outer corridors defined along the mean experimental curve (Gehre et al. 2009). The corridor score uses the ISO default scheme. The phase, magnitude, and slope scores are calculated using the methods of Enhanced Error Assessment of Time Histories (EEARTH) (Zhan et al. 2011). The phase score is calculated by identifying the maximum cross-correlation between the model and reference curve after time-shifting the model curve by discrete intervals. The slope and magnitude scores are then calculated from the time-shifted curve corresponding to maximum cross-correlation. The slope score is obtained by comparing the average slope in 1-ms intervals between the reference and model curves. For the magnitude score, dynamic time warping (DTW) is first applied to minimize the effects of phase on the magnitude difference (Barbat et al. 2013). Following the application of DTW, the one-norm of the difference between the curves is calculated for each time point and normalized by the one-norm of the shifted and warped reference curve: (1) $$\begin{array}{ccc}\hfill \mathrm{Total}\mathrm{ISO}& =& 0.6\times \left(\mathrm{sub}\mathrm{metric}\mathrm{average}\right)\hfill \\ & & +0.4\times \left(\mathrm{corridor}\mathrm{score}\right)\hfill \end{array}$$(1)

Results

Results of the input kinematics showed highly repeatable experiments. For the simulations, the means of the traces in Figure A10 (see online supplement) were used. Kinematic data from each of the models were compared against NBDL volunteer response corridors for each simulation (Figures A4–A6, see online supplement). The corridors for visualization are defined as the average volunteer response plus and minus the standard deviation.

The dotted black line represents the average volunteer response and the solid black lines are the corridors of one standard deviation. Angular accelerations are given in radians per second squared and linear accelerations were converted to g's to match previous publications. Due to the directionally controlled nature of the testing setups, only rotational accelerations in the key planes of motion were measured from the models. The rotational Y accelerations shown have been sign-adjusted for the –Y to correspond to the NBDL output coordinate system. Positive rotational Y acceleration corresponds to frontal rotation and positive rotational X acceleration corresponds to lateral rotation to the right.

All three models are shown to behave with a similar response for each of the simulations. The curves generally follow the shape of the volunteer data with some noted phasing of the peaks. Timing of the results is likely a function of the unknowns regarding the experimental setups, particularly related to the belting. The ISO scores for the comparisons of each simulation to the corresponding volunteer test are shown in Figure 1. All submetrics of the ISO score and corridor for each model simulation are reported in Table A3 (see online supplement). The total scores are reported in Figure 1.

Figure 1. Graphical comparison of ISO scores from the M50-O, M50-OS, and M50-OS+B in 15g frontal, 8g frontal, and 7g lateral simulation.

Discussion

The models capture the trends generally for head motion during the three simulations. In all models, the linear acceleration of the head, like the volunteer data, is bimodal. The precise time to peak and length of trough shows some variation against the volunteer data, but this may be a function of a number of factors: volunteer head mass, neck length, muscle strength and activation state, the belt system, and interaction with the buck. Motion of the upper thoracic vertebrae, specifically T1, has a clear impact on the kinematics of the head. It appears that T1 linear acceleration was slightly underestimated for all models in each of the three simulations.

The M50-OS and M50-OS+B behaved similarly in all simulations as predicted and corresponded fairly well to the M50-O. The M50-OS+B and M50-OS models lag in terms of lateral motion but otherwise performed well when compared to the M50-O. The results indicate for frontal studies in particular that the M50-OS+B could finish nearly nine simulations on the given hardware in the same time that it would take the M50-O to finish one, thus enabling further parametric study.

In taking the approach not to prescribe T1 motion, but rather to try to model the motion of the volunteers using restraints, there were significant details that were not specified from the study and could not be used to validate the belt system such as belt loads and seat reaction loads. Therefore, aspects of the belting system such as pretensioner values, pretensioner times, and belt properties were estimated to better match the kinematics of the volunteers. The sensitivity of these selections was evaluated. Variation of the shoulder belt locking time yielded the largest range in ISO scores of 13 percentage points, yet other belt properties showed a more modest effect. Muscle activation sensitivity found that earlier activation improved the overall ISO score, which is consistent with likely “braced” instructions for the subject. Kinematics were found to be least sensitive to variations in shoulder pretensioning time and scaling the seatbelt loading/unloading curves. Locking of the shoulder belt immediately after pretensioning yielded the lowest ISO comparison for angular head acceleration, suggesting that better coupling with the seat would improve that measure. Generally, accuracy of T1 acceleration and angular head acceleration improved as locking time increased to 100 ms, while the accuracy of linear head acceleration decreased as locking time increased.

Subjective grading scales have been suggested for analyzing ISO scores with R > 0.94 as excellent, 0.80 < R ≤ 0.94 as good, 0.58 < R ≤ 0.80 as fair, and R ≤ 0.58 as poor. Note that there is no universal grading scale, as applications for comparison may vary in difficulty. Fair and generally comparable ISO scores were found across models. Previous studies of whole-body kinematics and kinetics (ATD or human model), where the experimental setup is more meticulously documented, produced quantitative comparison scores of roughly 0.6 to 0.7 (Pietsch et al. 2016; Vavalle et al. 2015), and for component tests that can be greater (Davis et al. 2016). Thus, we consider these results to indicate reasonable biofidelity. It should be noted that Pietch et al. used CORA rather than ISO, which is another objective evaluation method.

A comparison of the ISO scores shows that generally each model earned a “good” rating for phase, “fair” rating for magnitude, “poor” rating for slope, and a “fair” to “good” rating for corridor score. It appears that the slope submetric was the primary factor in lower ISO scores for the three testing simulations. ISO has been observed to appropriately grade slope relation between curves with low noise. However, since the slope score is calculated by comparing the average slope of each signal in 1-ms intervals, noise in the model curve detrimentally affects the slope score (Barbat et al. 2013). Therefore, noisy acceleration signals that appear to track the shape of the reference curve earned low scores even though intuitively based observations show the shapes to be similar. This study did not explore filtering techniques to combat this issue. The bimodal nature of the curves also increases the difficulty of obtaining high scores for the three ISO parameters. While the slope scores were the lowest, it should be noted in all cases that the models matched the head kinematic trends. In agreement with previous studies (Decker et al. 2016), it appears that the simplified models correlate best with frontal simulations.

Some differences of the two approaches to model the neck in these models bear some discussion in this study. The M50-O and M50-OS+B neck systems contain the same bony geometry but have different implementations. The M50-O cervical spine is deformable whereas the M50-OS+B cervical spine is a series of rigid bodies. The M50-O neck musculature is comprised of 1D active Hill type muscle elements embedded within three-dimensional (3D) passive neck muscles, which, in turn, are enveloped in a soft-tissue layer. The M50-OS+B only has 1D active Hill type musculature and a passive soft-tissue envelope. Intervertebral discs are physically modeled in the M50-O and this is implemented with 1D beam ligaments with progressive failure. The simplified model, on the other, uses a system of springs and dampers between cervical vertebrae (Dibb 2011). Despite these differences, head and T1 kinematics are similar between the two models, suggesting various paths to capture the fundamental biomechanics. The biggest difference is seen in lateral loading, which we attribute to the relatively large elements of the soft-tissue envelope in the simplified models. Furthermore, the work from Dibb (2011) that the simplified neck is based on did not include kinematic validation for lateral impact.

There are a number of limitations of this study. Those related to the belt system have been discussed already. Furthermore, how the active musculature of the neck was defined was held the same as previous work, in lieu of volunteer data from this work for the braced condition. The sensitivity study indicated that earlier modeling of bracing improves the ISO scores. The range of values for initial head CG–T1 distance between the volunteers was not reported in the publications and thus attempts were made to match a mean (which was determined through digitization of a published plot), yet it was not possible to match perfectly given anthropometric differences of the M50 models.

The M50-OS+B provides researchers with the option of obtaining representative brain response in one-ninth of the time of the full detailed model. This approach will be used in future parametric studies. Furthermore, gains could be realized, given that the simplified model is more likely to provide reasonable run times on smaller computational systems that may be more widely available. With additional tuning of the neck system in the simplified model, particularly in lateral loading, the M50-OS+B can expand on its potential.

Funding

Work was supported by the Global Human Body Models Consortium, LLC, and NHSTA under GHBMC project WFU-005.

Acknowledgments

The authors gratefully acknowledge the contributions Body Region Centers of Excellence (COE) for development and regional validation of the M50-O model. The GHBMC BRM COEs are located at Wayne State University (Head COE, PI Liying Zhang), the University of Waterloo (Neck COE, PI Duane Cronin), the University of Virginia (Thorax, Pelvis, and Lower Extremity COE, co-PIs Matt Panzer, Rich Kent, and Jeff Crandall), and IFSTARR (Abdomen COE, PI Phillipe Beillas). All simulations were run on the DEAC cluster at Wake Forest University with support by Damien Valides and Adam Carlson. Scott Gayzik is a member of Elemance, LLC, which distributes academic and commercial licenses for the use of GHBMC-owned computational human body models.