Comparison of control strategies for the cervical muscles of an average female head-neck finite element model

Abstract

Objective: ViVA OpenHBM is the first open source Human Body Model (HBM) for crash safety assessment. It represents an average size (50th percentile) female and was created to assess whiplash protection systems in a car. To increase the biofidelity of the current model, further enhancements are being made by implementing muscle reflex response capabilities as cervical muscles alter the head and neck kinematics of the occupant during low-speed rear crashes. The objective of this study was to assess how different neck muscle activation control strategies affect head-neck kinematics in low speed rear impacts.

Methods: The VIVA OpenHBM head-neck model, previously validated to PMHS data, was used for this study. To represent the 34 cervical muscles, 129 beam elements with Hill-type material models were used. Two different muscle activation control strategies were implemented: a control strategy to mimic neural feedback from the vestibular system and a control strategy to represent displacement feedback from muscle spindles. To identify control gain values for these controller strategies, parameter calibrations were conducted using optimization. The objective of these optimizations was to match the head linear and angular displacements measured in volunteer tests.

Results: Muscle activation changed the head kinematics by reducing the peak linear displacements, as compared to the model without muscle activation. For the muscle activation model mimicking the human vestibular system, a good agreement was observed for the horizontal head translation. However, in the vertical direction there was a discrepancy of head kinematic response caused by buckling of the cervical spine. In the model with a control strategy that represents muscle spindle feedback, improvements in translational head kinematics were observed and less cervical spine buckling was observed. Although, the overall kinematic responses were better in the first strategy.

Conclusions: Both muscle control strategies improved the head kinematics compared to the passive model and comparable to the volunteer kinematics responses with overall better agreement achieved by the model with active muscles mimicking the human vestibular system.

Keywords:

• Finite element analysis
• active muscle strategies
• rear impact
• whiplash
• female

Introduction

Human Body Models (HBMs) are usually developed with a specific human anthropometry, usually representing 5, 50, or 95th percentile of the adult male population (Pheasant and Haslegrave 2006; Östh, Mendoza-Vazquez, Linder et al. 2017). The ViVA OpenHBM was recently developed by Östh, Mendoza-Vazquez, Linder et al. (2017) to represent a 50th percentile female intended for assessing whiplash protection systems in a car. The ViVA OpenHBM was developed to address a gap in gender and size representation of the HBMs models. Epidemiological data shows that females are more likely to experience whiplash injuries (Carlsson et. al 2014). HBMs are one resource for investigating the injury mechanism if they can represent the biomechanics of the event.

Cervical muscles can play a significant role in whiplash injury causation by altering the head kinematics as well as their interactions with other anatomical sites in the neck such as facet capsules (Siegmund et al. 2009). To further study the effects of muscle contractions on whiplash injuries, activation time and force time histories for the cervical muscles should be considered (Siegmund et al. 2002). Data from the literature as summarized by Mang et al. (2012) have shown that the activation of cervical muscles occurs about 50-100ms after the onset of vehicle acceleration. This muscle response is considered to occur early enough in a rear impact event to alter the head and neck kinematics and affect the risk of injury.

To evaluate the influence of cervical muscle responses on whiplash injuries induced by rear impacts, active cervical musculature must be included in the ViVA OpenHBM. Two strategies have been previously proposed and are summarized by Östh et al. (2015) namely open-loop and closed-loop control strategies. In the open-loop strategy, previous simulations, experimental tests, or optimization results are used to predefine time history of a model’s muscle activation. In the closed-loop strategy, a control loop feedback function governs the activation of muscles. The intent for either approach is to mimic the feedback mechanism of the human’s Central Nervous System (CNS).

Muscle control depends on voluntary motion task, feedback through the body’s internal and external sensory systems, and reflexive motions. In this study, two main sensory feedback systems have been identified as relevant for modeling muscle response in a rear impact. The vestibular system and muscle spindles give feedback to the CNS for stabilizing the human head-neck complex (Keshner 2009). The sensory information from the vestibular system is used to balance the head and body in space by sensing the rotational and translational motion of head (Keshner 2009). Muscle spindles, which are found in high concentrations in the deep neck muscles (Keshner 2009), provide information to the CNS for the head-on-trunk orientation by sensing muscle length and the changes in the muscle length (Keshner 2009). Both these systems are relevant for the loading condition in a rear impact where the torso is disturbed and the resulting head and neck motions are assumed to be minimized by reflexive mechanisms steered by these two systems.

Recently, several researchers have included active muscles with feedback control mimicking the human vestibular system (Nemirovsky and Rooij 2010; Meijer et al. 2012; Östh et al. 2012; Iwamoto and Nakahira 2015), muscle spindle feedback (Feller et al. 2016), or both vestibular and muscle spindles feedback (Happee et al. 2017; Ólafsdóttir 2017). Happee et al. (2017), validated their model against small periodic perturbations in the sagittal plane while Ólafsdóttir (2017) used isometric and gravitational loading in five directions. To the authors’ knowledge, there is no study implementing active muscle contraction with feedback control in rear impact scenarios. The contribution of different controller strategies remains unclear during crash situations, particularly for car occupants that are unaware of the impending impact. Understanding the contribution of each specific controller is important to understand to successfully mimic head-neck kinematics with high biofidelity.

The objective of this study was to assess how different neck muscle activation control strategies affect head-neck kinematics in low speed rear impacts on the 50^th percentile female human body model.

Methods

ViVA open human body model (ViVA OpenHBM) Head-Neck model

The previously validated open source head-neck model from VIVA OpenHBM (Östh Mendoza-Vazquez, Linder et al. 2017; Östh, Mendoza-Vazquez, Sato et al. 2017) was used in this study. The head of the ViVA head-neck model is represented by a rigid body with a mass equal to 3.58 kg following the EvaRID model head mass (Carlsson et al. 2014). The ligamentous cervical spine model consists of approximately 116000 elements and has been validated for quasi-static loading as well as in dynamic rear-impact simulations (Östh Mendoza-Vazquez, Linder et al. 2017; Östh, Mendoza-Vazquez, Sato et al. 2017). Cervical muscles in the ViVA Head-Neck model are modeled as one-dimensional Hill-type muscle elements in LS-Dyna. The muscle implementations are based on the properties from Borst et al. (2011). In total, 129 fascicles of muscles were added to the model to represent muscles in the neck. In addition, a simplified cervical spine model was also developed by Östh Mendoza-Vazquez, Linder et al. (2017) by replacing the intervertebral soft tissues with kinematics joints.

The head-neck model is run separate from the remaining anatomy to reduce computational demands. By using nodal constraints, the distal ends of soft tissue and the cervical muscle elements, normally connected to the upper torso, were fixed to move together with T1 of the model (Östh Mendoza-Vazquez, Linder et al. 2017). Objective rating result using CORA shows only small differences between the detailed and simplified neck models when implemented in whole body simulations (Östh Mendoza-Vazquez, Linder et al. 2017). In the present study, the ViVA OpenHBM head-neck with simplified cervical spine model was used. More details of the ViVA OpenHBM is available in Östh Mendoza-Vazquez, Linder et al. (2017), Östh, Mendoza-Vazquez, Sato et al. (2017).

Control strategy to mimic neural feedback from the vestibular system

A Proportional Integral Derivative (PID) closed loop control was used to mimic the neural feedback from the vestibular system. The vestibular system was assumed to be a simple control system that maintain the head’s horizontal orientation during the impact. This strategy was defined using the LS-Dyna PID control function (PIDCTL) (LSTC 2016) and adopted from Östh et al (2012) and Ólafsdóttir (2017). To mimic the human body’s vestibular system, the coordinates of two nodes (Head Center of Gravity node and T1 node) and a reference node fixed to the head were used to define the controller vector. A reference value for the vector was taken at the start of impact. The angle between the reference and controller vector during the impact time defines the error signal for the feedback system. This error signal was delayed using the LS-Dyna delay function (LSTC 2016) to mimic human transmission and processing time delays. The delayed error signals were input to the PID controller to compute a control signal. This signal was filtered using a muscle activation dynamics filter, scaled by a spatial tuning pattern (Ólafsdóttir et al. 2015), and then input to the corresponding muscle fibers to generate a contraction force request. This control strategy will be referred as Angular Positioned Feedback (APF) Control (Figure 1).

Figure 1. Implementation of Angular Positioned Feedback (APF) control and Muscle Length Feedback (MLF) Control.

Control strategy to mimic muscle spindle feedback

Muscle spindle feedback was assumed to stabilize head-on-trunk orientation using a feedback system that tries to keep the muscle pre-defined length constant. The closed loop control strategy adopted from Ólafsdóttir (2017), also based on LS-Dyna PID control function (LSTC 2016), was implemented. The initial muscle length was calculated from the insertion and the origin nodes for each muscle element at the start of impact. The muscle lengths were calculated at each sample time and compared to the initial length to generate an error signal. As in the APF controlled, the error signal was delayed to mimic neural processing delays. Again, a PID controller uses the error signal to generate a control signal which is passed through activation dynamics filters and a spatial tuning pattern (Ólafsdóttir et al. 2015) to generate muscle force requests. This control strategy will be referred as Muscle Length Feedback (MLF) control (Figure 1).

Simulation setup

PID controllers have different gains for the proportional (KP), derivative (KD), and integration (KI) elements of feedback system. Optimization simulations for each controller (APF & MLF) were conducted to identify specific gain values (KP, KI, KD), for each controller (Figure 2). These optimizations searched for PID gains that would allow the FE model to reproduce human response. Published volunteer data in low-speed rear impacts from Sato et al. (2014) were used as reference data. The volunteers were two 23 years old females, 164 cm and 162 cm tall, and weighting 50 kg and 46 kg respectively. The volunteers’ statures were quite close to the ViVA model (height 161.6 cm and weight 60.8 kg) although a bit lighter.

Figure 2. Simulation setup and sign convention.

The data was derived from Test Series 2 (Sato et al. 2014) which used a rigid seat with a 20-degree vertical seatback angle without any head restraint. The seat was mounted to a sled that was accelerated by a compressed spring. The velocity change was equal to 5.8 km/h with a peak acceleration of 42 m/s² (4.3 G).

In the current study, the time histories of average linear displacement (x- and z-directions) and angular displacement (y-direction) of the first thoracic vertebrae spine (T1) from two female volunteers were prescribed for the T1 motion of the model. The total duration of each simulation was 400 ms which includes 100 ms of dynamic settling for a gravitational acceleration equal to 9.81 m/s². The 100 ms of settling was required to allow the cervical spine to achieve static equilibrium in a gravity field and ensure that initial position of the model was as realistic as possible. It was assumed that this setup was comparable to a full body simulation.

The gains for both control strategies were calibrated by conducting an optimization-based parameter identification using LS-OPT version 5.2 (Stander et al. 2015). The average head Center of Gravity (C.G) linear and angular displacement of the volunteers were used as the objective of the optimization simulation. The optimization was limited to 10 iterations (each iteration consisted of seven simulations) due to calculation time requirements. The optimization simulation was also used to define the neural transmission and processing delays of each controller.

Software and computational environment

Pre- and Post-processing of the simulation models were conducted using LS-Prepost (version 4.5-x64) with the additional of OriginPro 2018b (64-bit) as post-processing software. The simulations were run using LS-Dyna version nine double precision (ls-dyna_mpp_d_r9_2_119543_x64_redhat54_ifort131_sse2_intelmpi-413 for the APF simulations and ls-dyna_mpp_d_BETA_R9_121624_centos65_intel131_intelmpi for the MLF simulations). The reason for using different a LS-Dyna binary for calibrating the MLF controller was because the MLF controller was very time consuming since each muscle element has a PID controller. The latter binary version allows for the calculation to be conducted every 100^th circle of time step instead of every time step. All simulation results were normalized after the settling phase of each simulation, therefore data analysis was started after 100 ms.

Quantitative evaluation

To quantify the kinematics response between the model and the volunteers, deterministic-based metrics evaluation was conducted using a software called RSVVP (Roadside Safety Verification and Validation Program) (Mongiardini et al. 2013). RSVVP software in the current study was used to evaluate the similarities in magnitude and similarities in shape of volunteer kinematics results and simulation kinematics results. Although, all sixteen available metrics in RSVVP software were used to evaluate the simulation results, only Sprague-Geers metric without any pre-processing steps (see Appendix 1 (see online supplement)) was analyzed since several metrics have quite similar formulations (Mongiardini et al. 2013). In addition, objective rating evaluation using Correlation Analysis (CORAplus) software 4.0.4 (PDB, Gaimersheim, Germany) was also conducted (Gehre et al. 2009). CORA default corridors of 5% for inner and 50% for outer limits were selected as corridor definitions. This method was combined with the CORA correlation method with the final score based on equal weight of those methods.

Results

APF and MLF control gains

Both optimization simulations were numerically stable and converged according to the predefined iterations number. The gains were successfully identified and shown in Table 1. For the APF controller the best proportional and derivative gains were equal to 0.60 (%contraction/rad) and 412.62 (%contraction/rad ms-1) and for the MLF controller the best gain values were equal to 0.45 (% contraction/mm) for the proportional gain and 6.71 (% contraction/mm ms-1) for the derivative gain. From previous experience, the integral gains were constrained to zero (Östh et al. 2012).

Table 1. PID Gains from calibration simulation.

Active muscle control strategy	Proportional gain KP	Integral gain KI	Derivative gain KD	Neural delay
Angular Positioned Feedback (APF)	0.60 (% contraction/rad)	0 (% contraction/rad ms)	412.62 (%contraction/rad ms^−1)	15.00 (ms)
Muscle Length Feedback (MLF)	0.45 (% contraction/mm)	0 (% contraction/mm ms)	6.71 (%contraction/mm ms^−1)	53.21 (ms)
Passive	All gains zero

Muscle activation

Figures 3 and Figure A1 (see online supplement). Appendix two present the head kinematics of the passive model along with models with APF and MLF muscle implementation. The reference time T = 0 s is after the model settling is complete and the crash loads are applied to T1. The muscle element time series data provided insight into the response times for the different controllers, such as the sternocleidomastoid muscle forces shown in Appendix 2 Figure A1(d) (see online supplement). The onset of muscle force with the APF control was first observed between 55 and 75ms while the MLF controller started muscle activation later, between 120 and 150ms. Muscle activation changed the head kinematics of the model in both translational x- and z-direction as well as in angular y-direction compared to the model without muscle activation (Figure 3).

Figure 3. Time-series comparison of model with APF, MLF and without muscle activation.

Head kinematics and displacement

The activation of the muscles in the controllers started to alter the head kinematics in x-direction from around 70 ms (Appendix 2 Figure A1.a (see online supplement)). A model with the APF controller closely followed volunteer’s head motion until around 250 ms before starting to deviate. The MLF controller started to change the model head kinematics from 150 ms after impact. Both controllers improved the horizontal head motion compared to the passive model (no muscle activation) although the APF controller gave the best agreement.

For vertical head motions (z-displacement), APF and MLF controllers reduced the excessive peak head motions observed for the passive model (Appendix 2 Figure A1.b (see online supplement)). No numerical model could capture the first minima of volunteer head motion at 100 ms. Note that the sign convention reports upward motion as negative. The models with active muscle control reduced the total vertical range of motion for the head but there were still discrepancies between the model and the volunteer responses.

Comparison of head y-angular displacements (Appendix 2 Figure A1.c (see online supplement)) shows only small differences between the active and passive modeling approaches until 100 ms. The passive model was able to follow the volunteer motion best until around 200 ms (Appendix 2 Figure A1.c (see online supplement)). Models with APF and MLF controllers started to constrain the head angular motion after 200 ms although they did not fully capture the volunteer response curves.

Metric evaluations

Evaluation metrics are presented for the three modeling approaches in Table 2. The error values in bold identify model with the smallest error for the Sprague-Geers metric error and the highest score for the CORA evaluation. The bold values identify the best performing model. For head x-displacement, the model with active muscle control greatly improved the agreement with volunteer data and the smallest comprehensive error (combination of magnitude and phase error) was obtained for the model with an APF controller. In head z-displacement, a similar trend was observed, although the error of the best result was around 57.8%. In head angular y-displacement, the smallest comprehensive error was obtained by the model with the APF controller, however when the model with a MLF controller was compared to the passive model, better agreement was observed in the passive model as indicated by smaller error value. The values based on CORA gave similar results, showing the model with APF controller had the best correlation with the volunteer responses.

Table 2. Sprague-Geers metrics evaluation and CORA rating evaluation of calibration simulation.

Rating Evaluation of Kinematics	Head linear x-displacement			Head linear z-displacement			Head rotation y-displacement
	Passive	APF	MLF	Passive	APF	MLF	Passive	APF	MLF
Sprague-Geers Magnitude (%)	28.5	−2.9	20.4	191.5	51.9	147.8	44.1	13.9	48.9
Sprague-Geers Phase (%)	1.3	3.3	3.3	18.4	25.5	21.3	3.5	8.8	4.4
Sprague-Geers Comprehensive (%)	28.6	4.4	20.6	192.4	57.8	149.3	44.3	16.4	49.1
CORA Evaluation	0.707	0.940	0.782	0.324	0.344	0.321	0.681	0.707	0.615

Discussion

Closed-loop PID controllers for controlling cervical muscles in a FE head-neck complex were successfully implemented in the VIVA OpenHBM. Optimization simulations could establish the controller gains that would best fit head displacement information of the model to volunteer data based on Sato et al. (2014).

In general, the muscle activations improved the agreement between volunteer and model kinematics by reducing the head’s peak linear displacements compared to the model without muscle activation. The head angular displacements were higher in the model with active cervical muscles compared to the passive model.

When the model kinematics were analyzed in detail, a good agreement with the volunteers’ head x-displacement was observed in the model with an APF controller. This was not unexpected since the main function of the APF controller was to maintain the head center of gravity (C.G) in the horizontal direction. However, as the neck muscles are recruited to keep the head steady in the x-direction, the vertical motion became worse due to cervical spine buckling caused by the increased muscle tension seen in Appendix 2 Figure A1.(e) (see online supplement) at time 120–130 ms when there are the highest intervertebral rotations.

The unsatisfactory agreement of z-kinematics and rotational y-kinematics in the model with APF controller could also be caused by the lower agreement of the passive model itself. In the development and validation of the passive ViVA Model, Östh et al. (2017) pointed out that the CORA evaluations for the cadaver kinematic and the full body ViVA model in y-rotation and z-translation were around 0.53 and 0.31, respectively, which is rated as unsatisfactory. These ratings were similar regardless if the simplified or detailed head-neck model was used. Less cervical spine buckling was observed in the model with the MLF controller although the head kinematics were not as well captured as in the model with APF control. Less cervical spine buckling in the model with MLF controller was a result of the controller trying to maintain the muscle length constant, without overshooting. The kinematic differences between active models and the volunteer response in head y-angular displacements were mainly caused by the sudden sliding that occurs between the cervical spine vertebrae when they rotate. The models used a contact definition between the cervical vertebrae without any friction.

The APF controller was observed to respond faster to changes of head orientation for rear impacts compared to the MLF controller. Therefore, the APF controller seems more critical in the current scenario as it doesn’t rely on the muscle length changes that might be slower to detect changes in head orientation. However, the MLF controller could play a central role for preventing neck buckling and ensuring physiological limits for muscle forces.

To improve the model agreement with volunteer data, it could be beneficial and more physiological to combine the APF and MLF controllers into one global controller as attempted by Happee et al. (2017) and Ólafsdóttir (2017). Happee et al. (2017) also found that the controller that mimics muscle spindle response is important to reduce the cervical spine buckling when using a vestibular based feedback sensor. Both of these studies used volunteer data from less violent loading and therefore less range of motion. The results from this study demonstrate that both APF and MLF controllers cannot fully capture the volunteer motion over the entire event separately. The combination of the APF and MLF controller is beyond the scope of the current study.

It was observed in the current study that the force activation times for relevant muscles in rear impact collisions were different for the APF and the MLF controllers. This was because the optimal neural processing time constants, obtained from the optimization simulations for APF and MLF, were different. A wide range of muscle onset delays was observed in various volunteer studies as summarized by Mang et al. (2012) and Siegmund et.al (2002) therefore, it might be worthy to investigate further regarding how it affects the head motion during rear impact collision. The APF controller exhibited the best performance based on the historic data.

The PID strategy for both the APF and MLF strategies used constant controller gains. This may not fully represent the physiological processes but the research is attempting to produce a robust control strategy that can be tuned with available data sources. The single gain value implies that the controllers are not time dependent and thus assumes that the CNS response at the start of the impact is the same as at the end. This approach also assumes all cervical muscles in the model (258 beam elements of cervical muscle in right and left side) use the same global gain value. This simplification the constrains each muscle to the same control signal. Different levels of muscle activation in the cervical spine were observed in several of volunteer studies (Ólafsdóttir et al. 2015; Siegmund et al. 2007). Although, this might be a limitation, it was addressed in the current study by the assigning spatial tuning pattern before the muscle signal goes to specific muscle activation card. The spatial tuning changes the magnitude of a muscles control signal based on its contribution to a specific motion. It could be beneficial to group the cervical muscles into several groups to identify several gain values instead of using only one global gain value for every muscle.

Limitations

In the current study, the average responses of two female volunteers seated in a rigid seat were used to define the gains of the present model. Since the responses were only based on two volunteer responses, it might not be enough for representing the kinematics of average female population. The VIVA OpenHBM is intended to reflect the 50^th percentile female and thus individual volunteer response curves were not the reference data of interest, instead comparison to a representative average of the population was used. The initial cervical spine alignment pattern in the model and the cervical spine kinematics were not thus analyzed.

The use of a rigid seat instead of automotive seat with head restraint might not replicate the real crash scenario. The compliance effect of the seat cushion on the muscle activation time needs to be established but was beyond the scope of this study. The use of a rigid, controlled seat made the modeling and verification process more focused and reduced confounding factors like seat back effects.

The use of global PID gains for all muscle elements was a pragmatic approach to reducing the complexity of the controller. Potential development of gain scheduling or more advances control system designs are under investigation.

Acknowledgments

The simulations were performed on resources at Chalmers Center for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC) and carried out at Vehicle and Traffic Safety Research Center at Chalmers (SAFER). Special thanks to Karin Brolin, Jonas Östh, Jona Marin Ólafsdóttir for valuable discussions at the beginning of this research. Also special thanks to Erik Brynskog for the discussion regarding LS-Dyna code implementation.

Additional information

Funding

This study was funded by the Swedish Governmental Agency for Innovation Systems (VINNOVA).