Characteristics of head frequency response in blunt impacts: a biomechanical modeling study

Existing evaluation criteria for head impact injuries are typically based on time-domain features, and less attention has been paid to head frequency responses for head impact injury assessment. The purpose of the current study is, therefore, to understand the characteristics of human body head frequency response in blunt impacts via finite element (FE) modeling and the wavelet packet analysis method. FE simulation results show that head frequency response in blunt impacts could be affected by the impact boundary condition. The head energy peak and its frequency increase with the increase in impact; a stiffer impact block is associated with a higher head energy peak, and a bigger impact block could result in a high proportion of the energy peak. Regression analysis indicates that only the head energy peak has a high correlation with exiting head injury criteria, which implies that the amplitude–frequency aggregation characteristic but not the frequency itself of the head acceleration response has predictability for head impact injury in blunt impacts. The findings of the current study may provide additional criteria for head impact injury evaluation and new ideas for head impact injury protection.


Introduction
Head impact injury is a major public health issue, and the prevention and treatment of such injuries have become a significant challenge in the medical field globally, as they have a severe impact on the human body, not only leading to varying degrees of behavioral and functional disorders but also having a high fatality rate (MacDonald et al., 2014;Popescu C et al., 2015).There are many sources of head impact injuries, including falls, traffic accidents, dropped objects from heights, and sports such as football and boxing, among which traffic accidents are the dominant source of head impact injuries (Organization, 2018).
There are two types of criteria commonly used for head impact injury evaluation: one is based on head kinematic and dynamic responses such as acceleration and velocity, and the other is based on brain tissue biomechanical responses.The former mainly includes the head injury criterion (HIC) (Administration, 1996), generalized acceleration model for brain injury threshold (GAMBIT) (Newman, 1986), rotational injury criterion (RIC) (Kimpara and Iwamoto, 2012), head impact power (HIP) (Newman and Shewchenko, 2000), and brain injury criterion (BrIC) (Takhounts et al., 2013).The latter mainly uses the cumulative strain damage measure (CSDM), intracranial pressure, maximum principal strain (MPS), and strain rate (Cloots et al., 2013;Wright and Ramesh, 2012).However, due to the existence of the resonance phenomenon, head impact injury may produce a response in some frequency bands with a certain vibration mode, resulting in the "amplification" of the deformation of brain injury (Fonville et al., 2022).For this case, the traditional head impact injury evaluation criteria based on the head response in the time domain make it difficult to comprehensively assess the mechanical response of the head under dynamic impact, and joint analysis of the response in the time and frequency domains is required to fully understand the characteristics of the head mechanical response signal.
The finite element (FE) modeling method, the reduced-order model method, and the skull-cerebrospinal fluid (CSF)-brain system fluid-solid coupling method are mainly used in studies of the frequency response of human body head mechanical response.El Baroudi and Razafimahery (2014), El Baroudi et al. (2012a), and El Baroudi et al. (2012b) developed a skull-cerebrospinal fluid-brain fluid-solid coupling model for modal analysis and found that the fundamental frequency of the head was 26.66 Hz.Laksari et al. (2015) proposed that the skull-brain dynamic can be approximated as an under-damped system and reported that the head exhibits significant resonant behavior at approximately 15 Hz ± 2.9 Hz (Laksari et al., 2015).Gabler et al. established a single-degree-offreedom mechanical model to represent the brain-skull system and found that the intrinsic frequency of the brain was 22.3 Hz-27.5 Hz, which was close to the pulse duration (36 ms-45 ms) of the brain resonance period, and the shape of the brain displacement depends on the magnitude of velocity and acceleration (Ghodrati Amiri and Asadi, 2009).Yang et al. (2017) found that the fundamental frequency of the head finite element model was approximately 35.25 Hz for different damping factors.Laksari et al. (2018) reconstructed a mild traumatic brain injury case caused by collisions in American football players through dynamic mode decomposition (DMD) and finite element analysis, obtaining a fundamental frequency of 28 Hz (Laksari et al., 2018).Gabler et al. (2018) established a second-order detuned brain reducedorder model based on three-degree-of-freedom coupling and demonstrated that the inherent frequency of the brain-skull system was in the range of 21.6 Hz-29.3Hz around the natural cycle of the brain.Recently, Fonville et al. (2022) extracted different free-end modes of the head FE model by eigenvalues and found that the fundamental frequency in the bound head mode was 22.3 Hz and the fundamental frequency in the free boundary brain tissue mode was 13.9 Hz (Fonville et al., 2022).However, most of the studies mentioned above focused on mild traumatic brain injury; the change in frequency-domain responses from mild to severe impact was not explored much.Furthermore, there is still a lack of combination analysis on frequency-domain response and traditional injury criteria.
The purpose of the current study is, therefore, to understand the characteristics of human body head frequency responses to blunt impacts and the correlation between head frequency responses and commonly used injury criteria.First, head blunt impacts under various boundary conditions were simulated using the FE modeling method.Then, the wavelet packet analysis method was employed to deal with the head acceleration signal for extracting frequencydomain characteristics and collecting energy.Finally, the characteristics of head frequency responses and their correlations with the commonly used kinematic-and biomechanical-based criteria were analyzed.
2 Materials and methods

FE simulation setup
The head model extracted from the THUMS AM50 V4.0 human body model (mentioned as the THUMS head model below) was used for head blunt impact simulation.The THUMS head model consists of 37,756 nodes and 49,598 elements; the components include the skull, dura mater, brain, cerebellum, and brainstem, where the brain is surrounded by a layer of CSF, and the inner cranial bone consists of a hard shell to represent the anatomical structure of the head and brain more accurately (Figure 1).This model has been validated for its biofidelity (Iwamoto et al., 2002;Iwamoto et al., 2015;Kimpara et al., 2006;Watanabe et al., 2012) and is widely used in head blunt injury analysis (Iwamoto et al., 2002;Iwamoto et al., 2015;Wang et al., 2022;Wang et al., 2023); particularly, the THUMS head model showed good agreement with cadaver impact test data in terms of response force, acceleration, relative displacement between the intracranial tissue and skull, brain pressure, and biomechanical response (Iwamoto et al., 2002;Iwamoto et al., 2015).
To simply simulate head blunt impacts under different boundary conditions, block-to-head impacts at various speeds were modeled using cylindrical impact blocks with different sizes and materials and the THUMS head model.A simulation matrix based on the full factorial test was defined using the parameters shown in Table 1.The detailed dimensions (min = 12.5 mm×12 mm, mid = 25 mm×12 mm, and max = 37.5 mm×12 mm) and impact locations (forehead center, right side near the wing point, and head top center) of the blocks are illustrated in Figure 2.
In the simulations, a fixed constraint was applied to the surface of the occipital foramen of the THUMS head model to limit the head motion; a surface-to-surface contact with dynamic and static friction coefficients of 0.3 was defined between the skin and impact block outer surface; and the direction of the impact speed was set in the normal direction of the block.Three impact speeds were selected from low to high based primarily on common vehicle-pedestrian and falling object crash speeds (Li et al., 2021).Three impact block sizes were selected based on the contact area with the head for studying the influences of different blunt objects and different contact areas on the frequencydomain response of the head in real collision accidents.Three materials (rubber, steel, and glass) were selected, mainly according to the vehicle parts that occupants and vulnerable road users may contact in accidents, e.g., the interior upholstery, A-pillar, and windshield, which were modeled using different material models and properties (Table 2) and can represent structures with different stiffnesses.Three impact positions were chosen to investigate the influences of different skull thicknesses and distances from different brain tissues on the frequency-domain response of the head, which can basically represent the common impact positions of the head given the symmetrical structure of the head.

Frequency response analysis
The wavelet packet transform (WPT) is one extension of the wavelet transform (WT) that provides complete level-by-level decomposition (Coifman and Wickerhauser, 1992).The WPT is commonly used to transform the sequence data on the time axis of non-stationary random signals (such as acceleration signals) into spectral data on time and frequency, which can provide information about the intensity of a non-stationary timedependent motion of interest at a specific frequency, and it is easier to obtain the energy changes of each frequency band (Gabler et al., 2019).Furthermore, wavelet packets have a high frequency resolution for nonlinear signals and can segment the frequency band to study the contribution of different frequency bands to impact energy, which has the ability to link the frequency band to head injury.Therefore, the WPT was used to study the frequency-domain characteristics of nonstationary random signals of the head generated by blunt impacts.
The principle of wavelet packets can be illustrated as follows.Assuming that the wavelet packet function is u n (t), it needs to satisfy the bi-scale equation where Z is a positive integer; h k and g k are two-scale functions, g k =(-1) k h(1-k), which have an orthogonal relationship.The sequence {u n (t)} constructed using Eqs 1, 2 is called the wavelet packet of the scaling function.When n = 0, u 0 (t) and u 1 (t) are the wavelet basis functions of the scale functions φ(t) and ψ(t), respectively.The wavelet packet decomposition program simply transforms the signal from time domain to frequency domain while maintaining equal energy (Lu et al., 2013).For an acceleration signal x(t) subjected to wavelet packet decomposition up to the ith level, the energy of each sub-band can be calculated using the following formula: and the formula for calculating the total energy is Appropriate wavelet basis is a key factor in the signal processing effect of wavelet packet analysis.Due to the fact that the skull-brain impact signal is a non-stationary random signal similar to seismic signals and mechanical fault signals, a previous study (Fonville et al., 2022) has shown that when the skull and brain are combined as a composite structure for analysis, the skull-brain frequency response characteristics are mainly characterized by low-frequency features and may exhibit certain abrupt signals.Therefore, an sym2 wavelet with a sampling frequency of 5000 Hz, which has tight support, poor regularity, low vanishing moment order, and orthogonality to maintain energy (Bianchi et al., 2015), was selected as the wavelet base in this paper.According to Shannon's sampling theorem (Shannon, 2001), the Nyquist frequency is 2500 Hz.The number of wavelet packet decomposition layers was set to 9, which means that the original signal is divided into 512 sub-bands in the entire frequency domain, and each sub-band represents a bandwidth of 4.88 Hz.

Data analysis
The acceleration signals output from the simulations were subjected to wavelet packet-based frequency response analysis, which was extracted from the contralateral skull at the center of  the head collision at a sampling frequency of 5000 Hz (i.e., sampling every 0.0002 s).This sampling frequency was defined by taking into account both the accuracy of the results and computing efficiency since a low sampling frequency may miss the key acceleration peaks and a high sampling rate may increase timing costs.Then, three potential head injury indicators in the frequency domain were calculated using Eqs 3, 4, which are (E ij ) max , (i,j) (Eij)max , and (E ij / E total ) max , denoting the energy of the sub-band with the maximum value over all sub-bands, the position of the sub-band with the maximum energy (multiplied by 4.88 Hz is the frequency band), and the proportion of (E ij ) max in the whole energy, respectively.These three frequency-domain parameters can effectively extract key frequency-domain features contained in the signal: (E ij ) max reflects the strength of the signal within that frequency band, (i,j) (Eij)max indicates the frequency band that contains the maximum energy in the signal, and (E ij /E total ) max denotes the concentration of frequency.Finally, the kinematic-based head injury criterion HIC and the biomechanical-based brain injury criterion MPS were also calculated for each simulation.Then, spectral analysis was conducted to understand the characteristics of head frequency responses; parametric analysis using the nonparametric test was performed to analyze the influences of impact boundary conditions on the outcome of (E ij ) max , (i,j) (Eij)max , (E ij / E total ) max , HIC, and MPS; and regression analysis based on linear fitting was carried out to study the correlation between frequencydomain head injury indicators and existing head/brain injury criteria HIC and MPS.

Spectral analysis
Figure 3 shows the typical energy-frequency response of head acceleration signals from different impact conditions, where the case of Front_Steel_Mid_6 m/s was regarded as the reference condition and each sub-figure shows the variation of a parameter in Table 1 based on this reference condition.Generally, the head energy is distributed in a wide frequency range of 0-2500 Hz but mainly within the frequency bands 0-500 Hz.Furthermore, four energy concentration frequency regions can be observed, i.e., 0-35 Hz, 60-350 Hz, 450-500 Hz, and 900-1200 Hz.It is also found that there is no obvious difference in the energy distribution of each frequency band when changing the impact location; the head energy at the low-frequency bands (<100 Hz) decreases with the increase in the stiffness of the impact block from rubber to glass to steel and also the impact velocity; on the contrary, an opposite trend was observed for the change in impact block size, where the head energy is more concentrated in the low-frequency bands (<100 Hz) when impacted by a larger block.
Figure 4 shows the distributions of brain MPS and skull von Mises stress under different impacts.Again, the case of Front_Steel_ Mid_6 m/s was regarded as the reference, and the variations of all parameters were represented in example cases (the caption below each sub-figure illustrates the impact boundaries).It could be found that the brain strain mainly occurs in the peripheral areas of brain Impact block size, impact position, and material.
Frontiers in Bioengineering and Biotechnology frontiersin.orgtissue, and the skull stress is mainly distributed in the impact block contact area and bone sutures.When changing the impact position from front to right to top, both MPS and von Mises stress gradually increase.The brain MPS and skull von Mises stress both increase with the increase in the stiffness, velocity, and size of the impact block.

Parametric analysis
Figure 5 shows the distribution of the normalized values for the time-domain and frequency-domain head injury indicators as a function of impact position, impact block material, impact velocity, and impact block size.For different impact positions, a wider fluctuation range and higher values of (i,j) (Eij)max and (E ij / E total ) max are observed in the top and front impact cases compared with the impacts at the right side, and the top impacts have the highest values of (E ij ) max , HIC, and MPS.When the material of the impact block shifts from rubber to glass to steel, both the fluctuation range and median value of the (Eij)max, HIC, and MPS increase, while the glass blocks lead to a generally lower (i,j) (Eij)max .A clear increasing trend was found for all indicators except (E ij /E total ) max when increasing the impact velocity.Generally, a smaller block results in a higher value and fluctuation range of the (i,j) (Eij)max and (E ij /E total ) max , while an opposite trend was found for the MPS.On the other hand, the middle-sized block induced a lower (E ij ) max and HIC.
Table 3 shows the results of the non-parametric test, which was performed to evaluate the statistical significance of the influences of the impact boundary conditions on the selected time-domain and frequency-domain head injury indicators.It is found that the impact position has no significant effect on all head injury indicators.However, the material of the impact block has a significant influence on the (E ij ) max , HIC, and MPS; the size of the impact block significantly affects (E ij /E total ) max and MPS values, while the impact velocity has a significant influence on all indicators except for (E ij /E total ) max .

Correlation analysis
Figure 6 shows the linear fitting results between the potential frequency-domain head injury indicators, (i,j) (Eij)max , (E ij /E total ) max , and (E ij ) max , and the HIC and MPS.It can be clearly seen that Typical energy-frequency response under different impact conditions.
Frontiers in Bioengineering and Biotechnology frontiersin.org(E ij ) max has a high correlation with the HIC (with a linear fitting R 2 value of 0.92) and a weaker linear relationship with the MPS (R 2 = 0.35), whereas (i,j) (Eij)max and (E ij /E total ) max show almost no linear correlation with both HIC and MPS.

Discussion
In this study, the wavelet packet processing program was employed for investigating the characteristics of head frequency response under blunt impacts.The analysis results (Figure 3) indicate that the impact energy observed on the head is mainly in the frequency bands 0-500 Hz; for few cases with a low impact speed (2 m/s) or a softer block (rubber), high head energy can be seen around the head fundamental frequency (14-35 Hz) reported in the literature (El Baroudi et al., 2012a;El Baroudi et al., 2012b;El Baroudi and Razafimahery, 2014;Fonville et al., 2022;Ghodrati Amiri and Asadi, 2009;Laksari et al., 2015;Laksari et al., 2018;Yang et al., 2017); however, for most cases, the head energy is concentrated in the frequency regions much higher than the fundamental frequency of the head.This finding might suggest that resonance issues could not be considered in the head blunt impacts.
The simulation results (Figure 5; Table 3) also show that the frequency-domain response of head energy under blunt impacts is significantly sensitive to the impact velocity and the material and size of the impact block, except for the impact location.Particularly, a higher impact velocity usually leads to a higher head energy peak (E ij ) max at a higher-frequency sub-band (i,j) (Eij)max , a stiffer impact block is associated with a higher head energy peak (E ij ) max , and a bigger impact block could result in a high proportion of the energy peak (over the total head energy) (E ij /E total ) max .It is easy to understand that higher impact velocity and a stiffer impact block (also higher density) could induce higher energy to the head, hence leading to a high energy peak (E ij ) max ; this could also be verified by observing the variation trend in HIC and MPS while changing the impact velocity and material of the impact block.However, the material characteristics of the impact block might have a different influence on the frequency compared to the amplitude, and hence, the frequency of the energy peak (i,j) (Eij)max is not sensitive to the stiffness of the impact block.On the other hand, the proportion of the energy peak over the total head energy (E ij /E total ) max is not affected by the impact boundary condition, given the fact that this parameter is affected by the value of both (E ij ) max and E total , which usually accompany an increase or decrease (a higher (E ij ) max is usually in case with a higher E total ).The impact position has no significant effect on head energy; this might be mainly because of the integrated structure of the skull, where different impact positions do not affect the absorption and propagation of impact energy.It is surprising that the size of the impact block has no significant effect on the head energy peak (E ij ) max .This may be due to the fact that as the size (and mass) of the impact block increases, the contact area between the head and the impact block also increases, resulting in more energy absorption and less pressure on the skull.However, a bigger impact block could induce more intense movements in the brain, resulting in a higher MPS.
Furthermore, the linear fitting results (Figure 6) between the frequency-domain parameters (i,j) (Eij)max , (E ij /E total ) max , and (E ij ) max and existing head injury criteria HIC and MPS indicate that only (E ij ) max has a linear correlation with HIC (R 2 = 0.92) and MPS (R 2 = 0.35).This may suggest that the head energy peak (E ij ) max has a certain contribution to head impact injury, particularly to skull fracture, as HIC was initially developed based on skull fracture risk (Coifman and Wickerhauser, 1992).The high linear correlation between (E ij ) max and HIC could also be easily understood by the fact that they are both based on the head acceleration signal.This finding may suggest that only the amplitude-frequency aggregation characteristic of the head acceleration response has head impact injury predictability under blunt impacts, while the frequency itself does not seem to be related to head injuries.
There are several limitations to the current work.First, the blunt impacts simulated in the current study were simplified as pure linear load; rotation may induce a significant influence, and future analysis may focus on real crash scenarios to further understand head frequency response.Second, the frequency response of THUMS has not been validated against biomechanical test data, although this model showed good injury-predictive capability in the literature.Finally, the defined frequency-domain injury predictors show a low correlation with the MPS, and other frequency-domain parameters would be focused on in further analysis.The current study is the first attempt to understand the characteristics of head frequency responses in blunt impacts of varying severity, where the FE human body models, wavelet packet frequency response analysis method, and statistical analysis method were employed.The simulation results show that the frequency-domain responses of head energy in blunt impacts could be affected by the impact boundary condition.The head energy peak and its frequency increase with the increase in impact; a stiffer impact block is associated with a higher head energy peak, and a bigger impact block could result in a high proportion of the energy peak.Regression analysis indicates that only the head energy peak has a high correlation with exiting head injury criteria, which implies that the amplitude-frequency aggregation characteristic but not the frequency itself of the head acceleration response has predictability for head impact injury in blunt impacts.The findings of the current study may provide additional criteria for head impact injury evaluation and new ideas for head impact injury protection.

FIGURE 4
FIGURE 4Typical distributions of brain MPS (left) and skull von Mises stress (right) under different impacts.

FIGURE 5
FIGURE 5Distribution of the normalized values for the time-domain and frequency-domain head injury indicators as a function of impact position, impact block material, impact velocity, and impact block size.

TABLE 1
Parameters for the definition of the simulation matrix.

TABLE 2
Parameters of the impact blocks in different materials.

TABLE 3
Non-parametric test results.
Item(i,j) (Eij)max (E ij /E total ) max (E ij ) max HIC MPS *A significant difference in the mean value at p < 0.05.Frontiers in Bioengineering and Biotechnology frontiersin.org