Abstract

1. Introduction

Traumatic brain injury (TBI) is a leading cause of morbidity and mortality in the USA that renders major socio-economic and public health burdens [1]. Indisputably, TBI is initiated through mechanical insults to the brain. Therefore, understanding the biomechanical mechanisms of TBI is critical to establish injury diagnostic criteria as well as to develop protective gear and rules to prevent and reduce the incidence, and severity of the injury.

Both linear (a_lin) and rotational (a_rot) accelerations contribute to head impact kinematics in real world. Therefore, both have been postulated as major TBI mechanisms. Elevated pressures resulting from a_lin-dominated translational/direct contact loading are typically related to focal brain injuries, which are also linked to skull fracture, epidural haematoma and contusion in moderate to severe brain injuries [2,3]. Early experimental studies reported that a higher peak pressure magnitude combined with a shorter duration or vice versa could both produce a concussive effect [4,5]. Subsequent studies on the significance of a_lin on cadavers [6] led to the use of the Wayne State Tolerance Curve (WSTC), Severity Index (SI [7]) and Head Injury Criterion (HIC [8]) to assess the risk and severity of head injury, including skull fracture and severe brain injury. Notably, HIC based on the resultant a_lin temporal profile is deemed as the official head injury standard, despite its well-recognized limitations for neither including a_rot nor considering variation in brain anatomy or impact direction [9]. Because these injury criteria are derived from a single resultant a_lin, which describes a point-wise rigid-body motion, they are unable to incorporate effects from non-rigid skull deformation that could also be important for brain pressures [10,11].

On the other hand, a_rot is considered to produce both focal and diffusive brain injuries [12]. The significance of a_rot was first hypothesized by Holbourn [13] and was later confirmed by experiments conducted on primates suggesting a_rot-induced strains responsible for concussive injury, diffuse axonal injury and subdural haematoma [14]. Because HIC does not include a_rot, its use in the mild/moderate spectrum of TBI has been criticized. Consequently, more recent injury metrics explicitly include a_rot or its derivatives for injury risk assessment. These include the generalized acceleration model for brain injury threshold (GAMBIT [15]), head impact power (HIP [16]) and the HIT severity profile (HITsp [17]). Analogous to HIC, injury metrics solely comprising a_rot have also been devised (e.g. the rotational injury criterion (RIC), power rotational head injury criterion (PRHIC) [18] and the brain injury criterion (BrIC [19])).

Unfortunately, these global kinematic variables and their variants do not directly relate actual brain injuries to tissue-level mechanical responses believed to initiate the injury [12,20,21]. In large part, therefore, no consensus has been reached on an appropriate kinematic injury metric to date. Further, there is still theoretical debate on the role of a_lin and its induced pressures on the risk of brain injury. Notably, injury metrics such as RIC, PRHIC [18] and BrIC [19] only consider a_rot but not a_lin. However, a recent study based on real-world on-field head impacts using a_rot and a_lin peak magnitudes suggests a_lin significantly improves a concussion risk function [22]. These controversial findings suggest a practical challenge in integrating theory and real-world experiment. Largely, this may be a result of using global head impact kinematic variables and their variants, as opposed to tissue-level brain mechanical responses [12,20,21], to assess the risk of brain injury.

Finite-element (FE) models of the human head are effective tools to estimate tissue-level responses using impact kinematics as input. However, current state-of-the-art computational schemes require a substantial computational cost (runtime and hardware) to simulate even a single impact. Therefore, they are not practical for real-world applications. Developing a computationally efficient, yet sophisticated numerical model with high predictive power remains a major challenge [23]. This challenge may be addressed using a pre-computation strategy. Using a pre-computed brain response atlas (pcBRA), (near) real-time estimation of brain strains was possible without significant loss of accuracy in the context of contact sports because of their similar temporal characteristics [24,25].

The goal of this study is to establish an analogous pcBRA to estimate a_lin-induced whole-brain pressures. Because a_lin-induced brain pressures are linearly related to a_lin magnitude and brain size [26], only the translational axis, Ω, (which is uniquely related to brain shape or anatomy) needs to be parametrized for a given head FE model, e.g. by the azimuth (θ) and elevation (α) angles. However, this requires 145 (13²–12 × 2) unique atlas solutions to sample one-half of the symmetric parametric space (assuming a 15° sampling density for θ and α, and further accounting for duplicated responses when α degenerates to ±90°, analogous to that in [24]). Here, we present an optimized and much more elegant approach to pre-compute only three atlas responses to achieve real-time pressure estimation, even on a laptop, for the entire impact duration (versus approx. 6 s for element-wise interpolation at one time point on a high-end multicore computer [24]) via linear scaling and superposition.

The pre-computed pressure and strain [24,25] response atlases, when combined, may help address the current debate on the roles of a_lin and a_rot on the risk of brain injury. This is an important step towards integrating theory and experiment in the field of brain injury prediction. Because of their real-time efficiency, they may enable future studies to focus on region-specific, tissue-level mechanical responses, rather than relying solely on global kinematic variables or their variants that have plagued early and contemporary TBI research to-date.

2. Material and methods

2.1. The Dartmouth head injury model

All brain responses were obtained using the Dartmouth head injury model (DHIM; figure 1). The DHIM is a subject-specific model created with high mesh quality and geometrical accuracy based on high-resolution magnetic resonance images of an athlete. Details of the model description, material properties and validation performances were reported previously [26–28]. Briefly, DHIM is composed of solid hexahedral and surface quadrilateral elements with a total of 101.4 k nodes and 115.2 k elements (56.6 k nodes and 55.1 k elements) and a combined mass of 4.562 kg (1.558 kg) for the whole head (brain). The average element size for the whole head and the brain is 3.2 ± 0.94 mm and 3.3 ± 0.79 mm, respectively. The DHIM uses a homogeneous Ogden hyperelastic material model with rate effects incorporated through linear viscoelasticity (‘average model’ in [29]) to characterize brain mechanical responses. The DHIM achieved an overall ‘good’ to ‘excellent’ validation against relative brain–skull displacement [30,31] and intracranial pressure responses [32,33] from cadaveric experiments, as well as full-field strain responses in a live human volunteer [34].

Open in a separate window

Figrue 1.

The DHIM showing colour-coded head exterior (a) and intracranial components (b), which also includes part of the spinal cord to improve model biofidelity in the inferior region. The x-, y- and z-axes of the model coordinate system corresponds to the posterior–anterior, right–left and inferior–superior direction, respectively. (Online version in colour.)

2.2. Pressure atlas responses

According to a recent dimensional analysis [26], brain pressure responses are essentially the result of contact forces generated from relative normal displacements of the brain at the brain–skull interface, which are uniquely determined by a_lin (magnitude and directionality), brain size and shape. Both a_lin and the resulting whole-brain pressure responses are functions of time. For simplicity, here we first consider a constant acceleration with a magnitude of a_lin and an arbitrary direction. The resulting whole-brain pressure, p, is hydrostatic. With a constant acceleration, the corresponding relative normal displacements of the brain at the brain–skull interface, d, is also a constant. It can be decomposed into three orthogonal components, d_x, d_y and d_z, respectively, along the three major axes of the head model coordinate system:

2.1

Because pressure, p, is a scalar value, we hypothesize that brain pressures corresponding to the three independent displacement components, p_x, p_y and p_z, respectively, can be simply superimposed to establish the brain pressure responses corresponding to d or effectively, the given a_lin:

2.2

Furthermore, as pressure is linearly related to the magnitude of a_lin and d, a set of baseline or atlas pressures corresponding to impacts along the three major axes, and respectively, can be pre-computed. These atlas pressures (referred to as ‘brain print’, analogous to ‘finger print’) can then be individually scaled linearly to obtain pressures corresponding to d_x, d_y and d_z, or effectively, p_x, p_y and p_z, respectively. Here, we have chosen a baseline peak a_lin magnitude of 100 g (for convenience and relevance to injury; 1 g = 9.8 m s⁻²) to establish the three baseline pressures using the DHIM (see below). This leads to the following set of equations

2.3

where a_x, a_y and a_z are the a_lin magnitude components along the three major axes (in g) at a given time. Combining equations (2.2) and (2.3), whole-brain pressure is readily obtained

2.4

It is important to note that the above ‘dimensional analysis’ considers a constant acceleration, or a_lin at a given instance in time. When the same analysis is applied to an arbitrary point in time for a real-world time-varying a_lin temporal profile, the resulting hydrostatic pressure essentially generates a first-order approximation, purposefully neglecting any potential dynamic inertial effect. Fortunately, the hydrostatic response (equation (2.4)) would be sufficiently accurate for most real-world impacts because brain pressure reaches an equilibrium almost instantly due to its high dilatational speed (e.g. with impulse duration greater than 2 ms according to [35]).

To establish the atlas hydrostatic pressures (i.e. and in equations (2.3) and (2.4)), a ‘ramp-and-hold’ a_lin profile (figure 2) was applied to the rigid skull along the three major axes. This allowed element-wise brain pressures within the 15–20 ms window to be averaged along the temporal direction to serve as the atlas responses. The choice of this particular temporal profile was to ensure a hydrostatic response and to avoid perturbations from the dynamic inertial effects (arrow in figure 2), as explained. The peak magnitude of a_lin (i.e. 100 g) was not important here, because the pressure responses would be linearly scaled during the stable ‘hold’ period (15–20 ms). Importantly, the resulting element-wise, hydrostatic pressures or ‘brain print’ served as atlas ‘modes’ for linear scaling and superposition (figure 3; equation (2.4)).

Open in a separate window

Figure 2.

A sample a_lin profile and the resulting P_coup and P_{c_coup} (pressure at the coup and contrecoup sites, respectively) for a frontal impact (i.e. along the x-axis). Element-wise pressure values within the 15–20 ms window (shaded area) were averaged along the temporal direction to serve as the atlas response. Pressure perturbations due to the dynamic inertial effects are evident (arrow). (Online version in colour.)

Open in a separate window

Figure 3.

(a–c) The three hydrostatic baseline or atlas pressures ( and referred to as ‘brain print’) corresponding to translational head impacts along the three major axes. These atlas pressures are linearly scaled and superimposed to directly assemble brain pressure responses at an arbitrary given time. (Online version in colour.)

3. Results

Figure 4 compares the estimated peak P_coup and P_{c_coup} with their directly simulated counterparts for all the cases combined. Except for case 3 (with an underestimation of 27.8% and 23.4% in peak P_coup and P_{c_coup}, respectively), the estimated and the directly simulated peak values nearly aligned. More quantitatively, table 2 summarizes the estimation performances in terms of the various metrics for P_coup (results for P_{c_coup} were similar and are not repeated). One representative case was chosen to illustrate the pressure estimation process (figure 5). Figure 6 graphically compares pressure distributions between the estimated and the directly simulated responses for three additional cases. For the remaining cases, the estimated P_coup temporal profiles are compared with their directly simulated counterparts, along with the resultant a_lin profiles (figure 7).

Open in a separate window

Figure 4.

Comparison between the estimated and the directly simulated peak P_coup and P_{c_coup} for all cases evaluated, demonstrating an excellent performance overall as most closely aligned with the diagonal line (except for one case with an underestimation of 27.8% and 23.4% in peak P_coup and P_{c_coup}, respectively; arrow). (Online version in colour.)

Open in a separate window

Figure 5.

Illustration of using three baseline or atlas responses ( and a–c) to construct a pressure distribution map for an arbitrary given time in a selected a_lin profile (case 5). The atlas responses are first linearly scaled according to equation (2.3) using the three a_lin magnitude components along the three major axes (d–f) and then superimposed (g), resulting in a virtually identical response relative to the directly simulated counterpart (h). The estimated P_coup and P_{c_coup} temporal profiles are compared with their directly simulated counterparts (i; time interval corresponding to HIC shaded). To improve visualization, only the resultant a_lin profile is shown (versus 3 d.f. used in pressure simulation and estimation). (Online version in colour.)

Open in a separate window

Figure 6.

Comparison between the estimated (a,d,g) and the directly simulated (b,e,h) whole-brain pressure distributions when P_coup reached its peak for three selected cases. The estimated pressure response temporal profile is nearly identical to the directly simulated counterpart, and is compared with that of the resultant a_lin (c,f,i). (Online version in colour.)

Open in a separate window

Figure 7.

Comparison between estimated pressure temporal profiles and their directly simulated counterparts for the remaining cases, along with the resultant a_lin profiles. (Online version in colour.)

4. Discussion

Using tissue-level brain responses may help resolve the current debate on the roles of a_lin and a_rot on the risk of brain injury. As tissue mechanical responses are believed to initiate the injury [12,20,21], they may enable a more reliable injury diagnosis than head impact kinematics alone. Obviously, this requires an efficient estimation of tissue responses in the first place for practical deployment, which is beyond the capability of current direct simulations using FE head models. In this study, we established a simple, yet effective technique to achieve real-time estimation—even on a laptop—of whole-brain, temporally resolved (i.e. whole a_lin profile) pressure responses in translational head impact. Essentially, this technique generates a ‘brain print’ along three orthogonal directions by pre-computing hydrostatic brain pressures corresponding to impacts along the three major axes to permit a direct superposition of whole-brain pressure responses. As only three baseline responses are required, individualized pressure response atlases can be easily generated without the need for mass or geometrical scaling in order to account for anatomical variations. The notion that an FE model is not needed when computing a_lin-induced brain pressures is not new (e.g. ‘hand calculation’ in [35]), albeit, baseline pressure responses from an FE model are still necessary. Because dynamic inertial effects are purposefully neglected, the pcBRA effectively generates a first-order approximation.

Nevertheless, the estimation scheme was sufficiently accurate for 11 out of 12 real-world 3 d.f. a_lin profiles according to a range of criteria (table 2). For the spatially local metrics measuring the accuracies of P_coup and P_{c_coup}, the average percentage errors for the peak values were within 4%. For their time-varying profiles, similar errors were found for the average nRMSE, either for the entire impact duration or the interval corresponding to HIC (i.e. temporally global or local). Their estimations were also considered ‘excellent’ based on correlation scores [43], regardless of the temporal window. The excellent performance was further confirmed by the spatially global ‘double-10%’ criterion that measured the accuracy of full-field pressure distributions (average V_0.1 of 1.9%, i.e. 1.9% of the whole brain experienced greater than or equal to 10% differences on average). Notably, these errors were far below those in real-world a_lin measurement (e.g. up to 17–31% according to [44,45]). We have also confirmed that the reference a_lin magnitude to pre-compute the atlas pressures was not important. When using 1 g instead of 100 g to pre-compute, the estimated peak P_coup (P_{c_coup}) pressure increased (decreased) by only 1.5% (0.6%) for case 5.

The failed case (case 3; table 2) had a much shorter impulse (approx. 1 ms) than the others (greater than or equal to 3–4 ms; figure 7c), although it also corresponded to the largest a_lin and HIC magnitudes. Significant underestimation was evident for both peak P_coup and P_{c_coup} (by 27.8% and 23.4%, respectively; albeit, still below the maximum a_lin measurement error; figure 4) or the entire brain (V_0.1 of 66%; table 2). Obviously, the underestimation was a result of the much more pronounced dynamic inertial effects activated by the short impulse (analogous to a mathematical ‘singularity’) not captured by the atlas responses. To verify this, we temporally scaled the 3 d.f. a_lin profile by five times while maintaining the same a_lin magnitude (leading to an impulse of approx. 5 ms with HIC of 33 093). Once again, the pcBRA achieved an excellent performance as the peak P_coup and P_{c_coup} percentage errors became 4.1% and 1.2% with V_0.1 of 0.0%. These results confirm previous reports that hydrostatic equilibrium predominates and wave effects are of no consequence in most real-world head impacts [46] (when impulse duration greater than 2 ms according to [35]). Additionally, the resulting peak P_coup and P_{c_coup} for the scaled profile (462.5 kPa and −489.4 kPa, respectively) far exceeded the thresholds for contusion (e.g. 173 kPa or 234 kPa corresponding to moderate or severe brain injury, i.e. AIS 3–4 or AIS 5–6, respectively [47]) or contrecoup cavitation (e.g. negative dilatational pressure of –100 kPa according to [48]). These findings suggest that the applicability of pcBRA does not necessarily depend on the magnitudes of a_lin or HIC, or the impact severity.

An important limitation of our study was the omission of brain dynamic inertial effects, which could be an important contributor to brain pressure responses with a very short a_lin impulse (e.g. less than 3–4 ms, figure 5i, figure 7a and c). While the majority of real-world impacts has an impulse duration greater than 2 ms [35] and the resulting brain pressures are essentially hydrostatic [46], it appears that a_lin-induced focal brain injuries such as subdural haematoma, epidural haematoma, subarachnoid haemorrhage and contusion often occur with relatively short duration (less than 5 ms) direct impacts [49]. For extremely short duration impacts (e.g. less than 1 ms, resulting from light projectile impacts), the brain pressures can increase significantly above the quasi-static estimates with large fluctuations at both the coup and contrecoup sites [50]. These observations suggest the importance to further compensate for dynamic inertial effects under these impact conditions in order to improve the estimation accuracy. This appears possible, though, e.g. by systematically probing the differences between the estimated and the directly simulated brain pressures subject to a number of pre-defined input variations (e.g. impulse duration and impact direction). The deterministic causal relationship could be explicitly modelled, making the compensation straightforward.

It should be recognized that the estimated pressures, even after compensating for the dynamic inertial effects, do not consider the potentially important non-rigid skull deformation, which was another limitation of our study. However, it should be noted that none of the existing injury metrics is able to incorporate effects from non-rigid skull deformation. Nevertheless, while non-rigid skull deformation does influence pressure responses [26,51–53], brain pressures are predominantly induced by the skull rigid-body acceleration when head protective gear (e.g. airbags and helmets) is used to extend impact duration, absorb impact energy and to minimize skull deformation [54]. On the other hand, when non-rigid skull deformation becomes more significant, a distributed impact force on the head or an impactor is necessary to explicitly model a deformable skull.

Interestingly, the pressure pcBRA derived from a rigid skull assumption may still serve as an upper bound for brain pressures when significant non-rigid skull deformation occurs. This is because the increased local contact area as a result of the skull deformation would lower P_coup (i.e. ‘over-estimation’; with an identical a_lin magnitude; albeit, uniquely determining a_lin magnitude may become challenging). The pressure over-estimation was parametrically investigated by iteratively adjusting the skull Young's modulus and impact force to maintain an identical a_lin magnitude in frontal impacts [26]. On the other hand, when significant non-rigid skull deformation is coupled with an extremely short impulse (e.g. case 3 in figure 7c), their respective significances (i.e. overestimate and underestimation, respectively) would compete. The resulting overall effects on pressure responses require further exploration.

Another limitation was that the pressure pcBRA only considers translational impacts but does not include a_rot-induced pressures, while both a_lin and a_rot contribute to head impact kinematics in the real world. This may explain the large differences in peak P_coup and P_{c_coup} compared with those found previously in which both a_lin and a_rot were applied, especially with large a_rot magnitudes (e.g. peak P_coup of 214.6 kPa reported here versus 138 kPa in [36] for case 5 with a_rot > 20 krad s⁻²). For smaller a_rot magnitudes (e.g. less than 5 krad s⁻² in helmeted contact sports), perturbations from a_rot-induced pressure appeared insignificant, and our study reported similar peak P_coup magnitudes (figures 6f and ​and7g;7g; 87.7 kPa and 73.9 kPa for cases 9 and 10 versus 81.5 kPa in [29]). These findings suggest that the pressure pcBRA may be directly applicable for translation-dominated impacts but a pressure compensation scheme may be necessary for impacts with large a_rot magnitudes. A systematic investigation of a_rot-induced pressures appears not to have been performed before, although a preliminary study suggests pressures from a_lin and a_rot can be superimposed [55]. However, the applicability of the pressure pcBRA when a significant a_rot component is present is probably debatable. This is because many consider a_rot or its resulting strain responses as the major injury mechanism in these cases [12], despite that combining both a_lin and a_rot peak magnitudes (presumably, the resulting pressure and strain responses) was found to improve the performance of an injury risk curve [22]. Therefore, it is conceivable that the pressure and strain pcBRAs be combined for optimal use in practice, especially for contact sports [24,25].

Certainly, these pressure response behaviours may be model dependent. Different CSF modelling strategies [56] and brain–skull boundary conditions [57] are known to significantly influence brain pressure responses. Therefore, it is unclear whether the real-time pressure estimation technique developed here is applicable to other head models (e.g. with an equation of state to simulate CSF [58]). However, brain pressure responses are insensitive to the brain nonlinear shear property as long as a sufficiently high bulk modulus is used to ensure a near incompressibility, which is well known [21] and was recently confirmed [26]. Regardless, the DHIM achieved an overall ‘good’ to ‘excellent’ performance when validating against cadaveric pressure responses, in addition to a similar performance validating against relative brain–skull displacements from cadaveric impacts and strain responses from a live human volunteer [26,28]. In comparison, many other models are typically not rated or validated against strain responses in a live human. Further, the DHIM also successfully confirmed an approximately 0.5 ms time delay of the pressure peak response relative to the a_lin peak magnitude (figure 7a) as experimentally observed [26,32], which, to the best of our knowledge, has not been reported before. These observations suggest high confidence in our findings.

In modern brain injury studies, theoretical debates still exist regarding the importance of a_lin and its resulting pressures on the risk of brain injury. Recent brain injury metrics such as RIC, PRHIC [18] and BrIC [19] are only composed of a_rot but do not include a_lin. However, on-field data based on real-world injuries using a_rot and a_lin peak magnitudes suggests a_lin does significantly improve the performance of a concussion risk function [22]. These controversies highlight the challenges when integrating theory and experiment in real-world applications. Our ‘brain print’ technique may serve as a useful and important tool to enable future brain injury studies to address this controversy using tissue-level brain responses as opposed to kinematic variables alone, even on a subject-specific basis. To allow further exploration of the use of the pre-computed ‘brain print’, we have made it freely available along with the associated analysis software for research use (see appendix A).

5. Conclusion

We have successfully established and validated a simple, yet effective technique based on three pre-computed atlas pressure responses for real-time, whole-brain, temporally resolved pressure estimation in translational head impact. The technique appeared sufficiently accurate and robust for a range of real-world impact/injury scenarios (with impulse duration greater than or equal to 3–4 ms). Therefore, it may allow future studies to focus on tissue-level brain responses, as opposed to global kinematic variables alone, to better understand the mechanism of brain injury.

Appendix A

The pressure pcBRA and a subset of strain pcBRA [24,25] were spatially resampled, and are freely available along with the associated analysis software for research use and evaluation, at http://www.thayer.dartmouth.edu/~songbai_ji/DHIM/.

References