Low-Rank Representation of Head Impact Kinematics: A Data-Driven Emulator

Head motion induced by impacts has been deemed as one of the most important measures in brain injury prediction, given that the vast majority of brain injury metrics use head kinematics as input. Recently, researchers have focused on using fast approaches, such as machine learning, to approximate brain deformation in real time for early brain injury diagnosis. However, training such models requires large number of kinematic measurements, and therefore data augmentation is required given the limited on-field measured data available. In this study we present a principal component analysis-based method that emulates an empirical low-rank substitution for head impact kinematics, while requiring low computational cost. In characterizing our existing data set of 537 head impacts, each consisting of 6 degrees of freedom measurements, we found that only a few modes, e.g., 15 in the case of angular velocity, is sufficient for accurate reconstruction of the entire data set. Furthermore, these modes are predominantly low frequency since over 70% of the angular velocity response can be captured by modes that have frequencies under 40 Hz. We compared our proposed method against existing impact parametrization methods and showed significantly better performance in injury prediction using a range of kinematic-based metrics—such as head injury criterion (HIC), rotational injury criterion (RIC), and brain injury metric (BrIC)—and brain tissue deformation-based metrics—such as brain angle metric (BAM), maximum principal strain (MPS), and axonal fiber strains (FS). In all cases, our approach reproduced injury metrics similar to the ground truth measurements with no significant difference, whereas the existing methods obtained significantly different (p < 0.01) values as well as substantial differences in injury classification sensitivity and specificity. This emulator will enable us to provide the necessary data augmentation to build a head impact kinematic data set of any size. The emulator and corresponding examples are available on our website¹.

Introduction

Traumatic brain injury (TBI) is one of the most debilitating health problems in our society today, with nearly two million new cases in the US every year (Taylor et al., 2017). The majority of these cases are considered mild, also known as concussion (Defense and Veterans Brain Injury Center, 2018). The substantial increase in reported concussions in contact sports (Selassie et al., 2013), together with the recent findings of increased long-term pathological changes (DeKosky et al., 2013), has sparked a public discussion and raised awareness about TBI. An important requirement is an accurate and objective diagnosis of concussions, which in turn could inform better protective equipment design and safer activities (Manoogian et al., 2006; Wu et al., 2014; Kuo et al., 2017; Kurt et al., 2017; Siegkas et al., 2019).

Head motion kinematics, including the rate, frequency, and direction of head's movement during collision, has been deemed as one of the most consequential metric in predicting brain injury. Historically, kinematic-based metrics such as head injury criterion (HIC) (NTSA, 1972), rotational injury criterion (RIC) (Kimpara and Iwamoto, 2012), and brain injury criterion (BrIC) (Takhounts et al., 2013) have been used to detect injury. These metrics are still widely used among researchers and are endorsed by safety regulating organizations such as the National Highway Traffic Safety Administration (NHTSA) (Laituri et al., 2016) and the National Operating Committee on Standards for Athletic Equipment (NOCSAE) (National Operating Committee on Standards for Athletic, 2012). More recently, brain tissue deformation-based metrics have been introduced that use head kinematics as input to computational models that can approximate the effect of head motion on brain displacement and deformation. These metrics either use simple discrete mechanical elements in lumped-parameter models, i.e., mass-spring-damper combinations, to give a rigid-body estimate of brain's relative motion with respect to the skull (Kornhauser, 1954; Low and Stalnaker, 1987; Laksari et al., 2015; Gabler et al., 2018b), or more complex finite element (FE) models with detailed geometry of the brain anatomy, which can simulate the local brain deformation and interaction with the stiff bony or membranous structures (Kleiven, 2013; Ji et al., 2014; Zhao et al., 2016). In the case of lumped models, brain angle metric (BAM), developed based on a data set of concussive and sub-concussive head impacts (Laksari et al., 2019), and in the case of FE models, maximum principal strain (MPS) and axonal fiber strain (FS) along the white matter axon fibers have been proposed as effective injury diagnosis metrics (Wu et al., 2019b).

Evidently, both for the kinematic-based and the brain deformation-based metrics, head impact kinematics play a major role. With the advent of wearable sensor technology, several groups have been collecting on-field head kinematic measurements during contact sports events (Hernandez et al., 2014; Laksari et al., 2018; Miller et al., 2018, 2019; Wu et al., 2018). However, despite these pioneering efforts, on-field head kinematic measurements are not widely available. As a result, researchers have resorted to simplifying parameterizations of head collisions as idealized biphasic acceleration impulses (Yoganandan et al., 2008; Ji and Zhao, 2014; Abderezaei et al., 2019). These biphasic impulses are commonly represented either by a triangle or half-sine and defined by two parameters: height and width constitute the magnitude and duration of a head impact impulse. The simplification of kinematic impulses serves the objective of emulating on-field kinematic data of a head impact with a few and manageable number of parameters to populate an otherwise infinite-dimensional loading space to investigate and establish a relation between head motion and brain injury. However, a potential disadvantage of these simplifications is overlooking valuable information that could prove detrimental in developing injury metrics. Therefore, it is paramount to understand the characteristics of real-world head impacts and whether we can accurately capture them through simplified approximations. Furthermore, advances in computational methods, including machine learning algorithms, have provided new and exciting avenues for fast and reliable prediction and diagnosis of brain injury. As a result, given the prohibitively high computational cost of current FE models, the biomechanics community has been trying to utilize such interpolative and machine learning techniques (Ji and Zhao, 2014; Cai et al., 2018; Wu et al., 2019a,b). However, a limitation of those techniques is the large number of kinematic data required to train these algorithms (in the order of thousands of head impacts; Wu et al., 2019a). Currently such a data set is not widely available. Thus, artificial augmentation of kinematic samples has been utilized as an alternative to satisfy that training data set requirements of such algorithms.

In this paper, we present a formal method to extract the most dominant features of on-field head impact kinematics from an existing data in the context of contact sports. We subsequently use the extracted features in order to construct an augmented data set that resembles the on-field measurements. We hypothesize that by using our method, based on principal component analysis (PCA), we will acquire more accurate injury predictions than current biphasic impulse approximations when compared against the ground truth measurements. Furthermore, we present a modal reconstruction technique that, despite using relatively few modes, can emulate a desired number of augmented head impacts that are statistically similar to the ground truth impact measurements.

Materials and Methods

In order to study the characteristics of head impact kinematics and the efficacy of simplified approximations, we used a previously-collected data set of 537 head impact kinematics measured during contact sports, including American football, boxing, and mixed martial arts (Hernandez et al., 2014; Laksari et al., 2018). For each impact, 6 degrees of freedom (DoF) kinematics—linear acceleration and angular velocity in the three anatomical directions—were collected at 1,000 Hz for 100 ms using a mouthguard instrumented with a triaxial accelerometer and a triaxial gyroscope (Hernandez et al., 2014). We construct three different reduced kinematics data sets to approximate the measured kinematics: (1) using principal component analysis (PCA), we decrease the dimensionality of the measured head impact kinematics to construct a low-rank kinematics data set, and (2–3) using previously proposed biphasic assumptions for acceleration impulses with acceleration magnitude and duration as the two variables, we construct biphasic data sets once for triangle (Tri) and once for half-sine (HS) approximations. We investigate the efficacy of each approximation by comparing its performance in detecting brain motion/deformation and injury prediction using three types of metrics: (1) kinematics-based injury metrics, including HIC, RIC, and BrIC, (2) brain angle metric (BAM), and (3) tissue deformation-based finite element injury metrics, including maximum principal strain (MPS) and axonal fiber strain (FS) in the whole brain (WB) and corpus callosum (CC). We compare the performance of each approximation against the ground truth (GT) data described above. We perform power spectral density (PSD) analysis on the derived temporal modes to obtain their predominant frequencies. These values are given by the maximum power spectral density (PSD) values of each mode. Finally, we present a modal reconstruction method for emulating augmented head impacts.

Dimension Reduction Through Principal Component Analysis

Our goal is to exploit the correlations between different measurements and find a reduced representation for three quantities of interest (QoIs), including: linear acceleration, angular velocity, and angular acceleration in each anatomical direction. In the case of linear kinematics, anterior-posterior, inferior-superior, and lateral directions are considered and in the case of angular kinematics, axial, coronal, and sagittal directions are considered as separate QoIs. To this end, we apply principal component analysis (PCA) to our data set. For each QoI, we form a data matrix X_m×n, where n = 537 is the number of measured head impacts and m = 100 is the number of time steps, and X = [x₁|x₂|…|x_n], where each column represents the measured QoI for a particular head impact and each row represents the time instance of the measurement. To perform PCA, we compute the singular value decomposition (SVD) of the data matrix: X = UΣY^T, where U = [u₁|u₂|…|u_n] are a set of orthonormal modes, i.e., $u_i^T{u}_j={\delta }_{ij}$, Σ = diag(σ₁, σ₂, …, σ_n) is a diagonal matrix, where σ₁ ≥ σ₂ ≥ ⋯ ≥ σ_n are the singular values, and Y = [y₁|y₂|…|y_n] are the uncorrelated linear components, i.e., $y_i^T{y}_j={\delta }_{ij}$ with the joint probability distribution function (PDF) of p(y₁, y₂, …, y_n). This gives an ordered array of the modal contributions inherent to head impacts response. Finally, we perform power spectral density (PSD) analysis on the derived modes u_i to obtain their predominant frequencies. These values are given by the maximum PSD values of each mode.

A reduced representation of the impact data is obtained by: $X\simeq \sum _{i=1}^k{\sigma }_i{u}_i{y}_i^T$. To quantify the performance of the reduction, we introduce:

$\begin{array}{ll}\eta (k)=\frac{({\displaystyle {\sum }_{i=1}^k}{\sigma }_i)}{({\displaystyle {\sum }_{i=1}^n}{\sigma }_i)},& (1)\end{array}$

Data-Driven Emulator

Once we extract the modal characteristics of the head impact data set, a reduced emulator of the head impact kinematics is obtained by truncating to k PCA modes:

$\begin{array}{ll}{x}^{*}={\displaystyle \sum _{i=1}^k}{\sigma }_i{y}_i^{*}{u}_i,& (2)\end{array}$

where $(y_1^{*},y_2^{*},\dots ,y_k^{*})$ is a random point with k components drawn from the marginal PDF of p(y₁, y₂, …, y_k). Equation (2) can be used as an emulator for producing new time series (x^*) for each of the QoIs that are nearly indistinguishable from the ground truth head impact kinematics measurements. Our approach can be interpreted as a stochastic dimension reduction technique and it is a special case of dimension reduction with time-dependent modes (Sapsis and Lermusiaux, 2009; Babaee and Sapsis, 2016; Babaee et al., 2017; Babaee, 2019; Patil and Babaee, 2020), in which the loading is a random input and the QoIs are random output. In order to ensure high quality augmented data, we establish features to statistically compare new time series with the GT, including parameters such as impact time, peak, and duration of impacts, and orthonormal projections of PCA modes. We will provide this emulator as a stand-alone program that allow users to build low-rank head kinematics data sets with various approximation levels (k) and emulate a pre-defined number of impacts.

Parameterizing Biphasic Impulse Profiles

In order to compare the performance of previously proposed biphasic models for angular and linear acceleration impulses (Ji and Zhao, 2014; Abderezaei et al., 2019) against our low-rank PCA approximation, we reconstruct biphasic triangle (tri) and half-sine (hs) representations of the 537 head impacts described above. First, the maximum absolute value of the angular acceleration profile (α_M)—computed by differentiating angular velocity measurements using a first-order forward divided difference method (two points)—was identified, including the time of peak (t_M). The impact duration (Δt) is defined as the time interval on either side of t_M. The boundary of this interval is defined as where the sign or the convexity of the acceleration profile (whichever comes first) changes. Convexity changes are computed based on the second derivative test using a common three-point stencil central finite difference derivative. This process is described in Figure 1, where the impact duration (Δt = t₁ − t₀) is the time elapsed between the initiation time (t₀ < t_M) and the completion time (t₁ > t_M) of impact. In the cases where t_M → 0 ms or t_M → 100 ms, since it is not possible to define t₀ or t₁ correctly, only half of the simplified pulse was created. Finally, the change in velocity (Δω) was computed as the area under the acceleration impulse ($\Delta \omega ={\int }_{t_0}^{t_1}\alpha (t)dt$), and the corresponding acceleration magnitudes for the triangle (α_tri) and half-sine (α_hs) approximations were calculated through:

$\begin{array}{ll}{\alpha }_{tri}=2\frac{\Delta \omega }{\Delta t},{\alpha }_{hs}=\frac{\pi \Delta \omega }{2\Delta t}.& (3)\end{array}$

Angular velocity pulses were then computed through direct temporal integration of low-rank angular acceleration pulses.

FIGURE 1

Figure 1. Constructing triangular and half-sine biphasic impulses (Ji and Zhao, 2014; Abderezaei et al., 2019): a representative angular acceleration trace in the coronal plane, magnified to better represent initiation time (t₀), completion time (t₁), time of peak (t_M), and peak acceleration (α_Tri and α_HS).

Accuracy of Injury Prediction Metrics

Kinematics-Based Injury Metrics

We used HIC₁₅, RIC₃₆, and BrIC to compare the performance of our proposed PCA reduction against the triangle and half-sine biphasic signals, i.e., triangle and half-sine approximations. To this end, we used previously published injury threshold values: (1) for HIC₁₅, values of 240 and 667 have been reported as 50% risk of concussion (Newman and Shewchenko, 2000) and skull fracture (Marjoux et al., 2008), respectively; (2) for RIC₃₆, a value of 10.3 × 10⁶ is reported as 50% risk of concussion (Kimpara and Iwamoto, 2012); and (3) for BrIC, a value of 0.5 constitutes a 50% concussion risk (Takhounts et al., 2013). We used these thresholds to assess the performance of each reduction approach in providing injury predictions in terms of sensitivity and specificity with respect to the ground truth measurements.

Injury thresholds were used as indicators, defining true positives (above the threshold) and negatives (below the threshold), while the predictive value of each impulse approximation was compared against the GT measurement.

Brain Angle Injury Metric

We further compared the performance of each approximation using injury criterion based on lumped-parameter models of the head. These models generally consider simplifying assumptions: skull and brain are considered rigid bodies and relative motion between the two represents a form of deformation and injury, and the compliance of the brain-skull interface such as the effect of bridging veins, dura and pia maters is represented by linear spring and damper elements (Kornhauser, 1954; Low and Stalnaker, 1987; Gabler et al., 2018a). Given the head kinematics as the base excitation input, these models can estimate the relative motion of brain and skull, particularly the angular motion since that has been seen as the more consequential type of motion (Sullivan et al., 2015). Recently, brain angle metric (BAM) was developed based on the characteristics of human brain and skull in finite element simulations, and validated against observed concussive and sub-concussive head impacts (Laksari et al., 2019). We compute BAM for each kinematic approximation (PCA, triangle, and half-sine).

Tissue Deformation-Based Injury Metrics

As a final step in studying the efficacy of the different kinematic approximations, we compared the performance of each approximation in predicting the tissue-level deformation metrics using finite element simulations, including maximum principal strain (MPS) in the whole brain (WB) and in the corpus callosum (CC) region, as well as axonal fiber strains (FS) in the corpus callosum region, which have all been proposed as predictive tissue-level metrics for injury classification (Kleiven, 2007; Ji et al., 2014; Laksari et al., 2018) (Figure 2). Recently a convolutional neural network (CNN) was developed based on pre-trained FE simulations based on the Worcester Head Injury Model (WHIM) (Zhao et al., 2017). This CNN method uses angular velocity data as input to approximate the regional brain deformations, i.e., maximum principal and axonal fiber strain (Wu et al., 2019a).

FIGURE 2

Figure 2. Representation of strains metrics with highlighted regions: (A) maximum principal strain (MPS) in the whole brain (WB), (B) maximum principal strain (MPS) in the corpus callosum (CC), (C) fiber strain (FS) in the corpus callosum (CC). Fiber colors represents directions: inferior–superior (blue), lateral (red), and anterior–posterior (green). Images were generated using 3DSlicer from ATLAS-based anatomical representation in FreeSurfer (Zhang et al., 2018).

Injury Metric Error Analysis

Having simulated the injury metrics for each head impact measurement, we calculated the corresponding injury metric (metric_GT) and the injury metric estimated by the kinematics approximation (metric_approx) using the equation below:

$\begin{array}{ll}error=\left|\frac{{\text{metric}}_{approx}-{\text{metric}}_{GT}}{{\text{metric}}_{GT}}\right|\times 100.& (4)\end{array}$

Subsequently, we performed Friedman test (Daniel, 1990) with a p-value of 0.01 (MATLAB, friedman) to show significant differences between each approximated impulse and the ground truth. We also performed sensitivity and specificity analysis to provide an estimate for the efficacy of approximating the metrics for injury diagnosis. For this purpose, we used previously published values for 50% risk of concussion, including MPS_WB = 0.2 (Patton et al., 2012), MPS_CC = 0.2 (Kleiven, 2007), and FS_CC = 0.074 (Giordano and Kleiven, 2014).

Results

Dimension Reduction Through PCA

We performed PCA on the measured kinematics data for the QoIs, i.e., linear acceleration, angular velocity and angular acceleration in each anatomical direction. We used the reduction criterion of η = 0.90, as defined in Equation (1), for all these cases. In the case of angular velocity, the minimum number of modes that satisfies this reduction criterion is k = 13, 15, and 10 modes for coronal, sagittal, and axial directions, respectively. In Figure 3 (top row), the PCA reconstruction of angular velocity in three anatomical directions for a sample case is shown. The sample case was chosen randomly from the 537 cases and it is represented by a column of the data matrix X. The ground truth measurement for the sample case as well as the reconstructed impulses with different levels of reduction are shown. It is clear that the 15-mode reduction yields a satisfactory reconstruction. In Figure 3 (bottom row), the individual and cumulative contribution of PCA modes are shown for the entire kinematics data set. These results demonstrate that with a relatively small number of PCA modes an accurate approximation of the head kinematic measurements can be achieved. As an additional analysis, in order to determine the convergence of this method, we performed PCA with several randomly selected subsets of the 537 ground truth measurements (with 100, 200, 300, 400, and 500 samples) and investigated the number of modes required to satisfy the η > 0.90 criterion. The results show that the minimum number of modes slightly grows with subsets size but levels off below 500 cases, indicating that our data set of 537 could be sufficient to reliably reconstruct a head impact data set (see Supplementary Figure 8).

FIGURE 3

Figure 3. (Top) Low-rank reconstruction of angular velocity using k = 3, 5, 15 PCA modes. (Bottom) Individual and cumulative contribution of PCA modes for angular velocity reconstruction. Columns from left to right show results for coronal, sagittal, and axial directions, respectively.

Data-Driven Emulator

The first three temporal modes show a classic modal behavior with 1, 2, and 3 peaks in all three anatomical directions (Figure 4). Together, these first three modes capture nearly half of the total angular velocity response, and with each additional mode, we can reconstruct a closer approximation with the ground truth. For more analysis of the modes (see Supplementary Figures 2, 4). To further study the distribution of these modal approximations, we show the orthonormal projection of the first five PCA modes (y₂, …, y₅) against the first and most energetic mode (y₁) of the PCA data (black circles and bars in Figure 5). It is clear that the modes follow a Gaussian distribution, which would be an important consideration for emulating more data points.

FIGURE 4

Figure 4. The three most energetic temporal modes for angular velocity for the entire 537 head impact data set.

FIGURE 5

Figure 5. Distribution graphs for second to fifth principal components (y₂, …, y₅) projected on to the first principal component (y₁) for angular velocity in the sagittal direction. PCA data (in black) follows a Gaussian distribution. Emulated data (red) is generated performing a Gaussian random number generator based on mean and variances from the PCA modes.

To show the performance of our data-driven emulator, we reproduced an additional k = 537 head impact cases by randomizing the $y_k^{*}$ columns each with the same mean and variance as the original Y matrix (Equation 2). Projection of the orthonormal vectors $y_i^{*}$ on $y_1^{*}$ for an augmented data set (red circles and bars in Figure 5) shows similar distribution as the ground truth data. In addition, features of the augmented data generated on our emulator, such as duration and peak distributions, have not significant differences with respect those of GT (see Supplementary Figures 5–7). Thus, our emulator is able to successfully generate reliable kinematics sets of data based on real on-field measurements. A copy of the head impact kinematics emulator will be available on our website (Arrué, 2020).

Natural Frequencies

The contribution of the most dominant frequencies, obtained through PSD criterion, are displayed in Figure 6. In general, low frequencies interval such from 10 to 40 Hz have the highest predominance for each parameter. Notably, the cumulative contribution for rotational velocity is progressively decreasing with increasing frequency.

FIGURE 6

Figure 6. Natural frequencies and their contribution to head motion kinematics.

Parameterizing Biphasic Impulse Profiles

Using the criteria described above, we fitted triangle and half-sine analog pulses to the ground truth kinematics measurements in order to parameterize the rotational acceleration magnitude and duration for each head impact. As a result, we derived 537 analog impulses in the three anatomical directions for both triangle and half-sine approximations. The results are presented in Table 1, where the rotational acceleration magnitude and duration average and standard errors of the mean are given for the ground truth and the impulses approximations.

TABLE 1

Table 1. Mean and standard error of acceleration magnitude and duration for ground truth and three approximations (PCA, triangle, half-sine).

Accuracy for Injury Prediction Metrics

Kinematic-Based Injury Metrics

The injury metrics HIC₁₅, RIC₃₆, and BrIC were computed for every model and the ground truth. Figure 7 shows the distribution of all samples and the corresponding concussion and skull fracture thresholds. The PCA predictions showed similar mean and standard deviations for HIC and RIC (38.20 ± 139.55 and 1.89 × 10⁶ ± 7.66 × 10⁶, respectively) as the ground truth (38.49 ± 140.13 and 2.26 × 10⁶ ± 9.04 × 10⁶); however, there is a significant difference between the ground truth predictions and the triangle (15.74 ± 41.71 and 3.98 × 10⁵ ± 1.06 × 10⁷) and half-sine (14.80 ± 39.31 and 3.86 × 10⁵ ± 1.02 × 10⁶) impulses.

FIGURE 7

Figure 7. Computed HIC₁₅, RIC₃₆, and BrIC for each of the ground truth (GT) and the three approximated data sets: PCA, triangle (Tri), and half-sine (HS). The circles represent each sample, the solid red line represents the mean, and the blue and red regions show the standard deviation and standard error, respectively. The solid blue and dashed red lines represent 50% risk of concussion and skull fracture, respectively. Significant differences are indicated with (*) for p < 0.01. HIC and RIC graphs are on log-scale.

Additionally, whereas the PCA-based impulses showed accurate predictions compared to the ground truth, we observed that the biphasic approximations either under-predicted injury (higher number of false negatives) in terms of HIC and RIC, or over-predicted injury (higher number of false positives) in terms of BrIC (Figure 7).

To better illustrate this, we performed sensitivity and specificity analysis for injury classification with respect to the ground truth, where the PCA-based signals showed high predictive performance compared to the biphasic impulses (Table 2).

TABLE 2

Table 2. Sensitivity and specificity of kinematics-based injury metrics for HIC₁₅, RIC₃₆, and BrIC for the three approximations: PCA-based method, triangles (Tri), and half-sine (HS) compared to the ground truth, considering thresholds of 50% risk of injuries.

Brain Angle Injury Metric

We computed 3DoF relative brain angles using the lumped model proposed in Laksari et al. (2019) to obtain the maximum resultant relative brain angle as a result of each head impact based on the ground truth kinematics and each of the three approximations. In coronal and axial directions, triangular and half-sine approximations gave significantly lower predictions for the brain angle metric, whereas the PCA modes showed no statistically significant difference from the ground truth (Figure 8). We observed substantially smaller approximation errors (Equation 4) for the PCA approach (~3%) compared to the two biphasic approximations (~25%) (Figure 8).

FIGURE 8

Figure 8. Computed brain angle metric (BAM) values for each impact case based on ground truth and the corresponding PCA and biphasic approximations. The significant differences are indicated with (*) for p < 0.01. The error comparison between models with means and SEM of each data set is shown on the right.

Tissue Deformation-Based Injury Metrics

Similar to brain angle metric, we calculated the errors between the ground truth strain metrics and those of low-rank approximations. As can be seen in Figure 9, the PCA approach closely follows the ground truth simulation results in all three strain metrics: maximum principal strain (MPS) in the whole brain (WB) and corpus callosum (CC), as well as axonal fiber strain (FS) in the corpus callosum. The strains computed for the biphasic impulses significantly differ from the ground truth values and successively provide higher errors as we include region-specific and morphological information. Furthermore, based on the 50% concussion risk thresholds mentioned above, the PCA approach provides substantially higher sensitivity and specificity for injury classification (Table 3).

FIGURE 9

Figure 9. Simulated strain metrics for ground truth head impact measurements and the corresponding PCA and biphasic approximations. The blue line shows the 50% risk of concussion for each metric. The significant differences are indicated with (*) for p < 0.01. Also the mean and standard error of the mean are shown on the right for strain estimation errors.

TABLE 3

Table 3. Sensitivity and specificity of strain metrics for maximum principal strain (MPS) for whole brain (WB) and corpus callosum (CC), as well as fiber strain (FS) for CC compared to the ground truth, with respect the threshold of 50% risk of concussion.

Discussion

In this study we provide a formal approach for reducing the dimensionality of head impact kinematics in contact sports settings. We first derived the most important modes contributing to the head kinematics through principal component analysis and then compared those to existing methods that approximate head kinematics with simple biphasic impulses. We show that the modal decomposition approach can capture the kinematic behavior of the head with better accuracy and provide better approximations of brain deformation and injury classification. This analysis confirms that although head kinematics during head collisions span a wide range of magnitudes and frequencies (Laksari et al., 2018), we can accurately capture the impact head kinematics by using only a relatively small number of modes. The low-rank database constructed based on PCA analysis require only 15 modes to build the ground truth angular velocity kinematics with over 90% accuracy as well as accurately capture the predictive value of head impact kinematics using a variety of injury metrics.

A major advantage of our approach is that with the acquired modes above, we are able to emulate a head impact data set with any given number of cases without needing access to the ground truth measurements. This emulated data set would closely replicate head impacts measured by on-field wearable sensors that constitute current state of the art.

The advantage of this low-rank emulator, in addition to its computational efficiency, is that it avoids simplifying assumptions for the shape of acceleration impulses and only uses empirical measurements. In contrast, the conventional biphasic assumption for head impacts as simple impulses with only two variables, i.e., acceleration magnitude and duration, falls short in providing accurate estimates. This effect is more pronounced for acceleration impulses that are more variable due to the time derivative but is true for velocity profiles as well. This apparent lack of accuracy in injury prediction in the biphasic approximations might be due to the fact that the biphasic triangle and half-sine signals are built using acceleration signals and then integrated to give the corresponding velocity profiles. Since there is no restitutive deceleration for these impulses, angular velocity eventually becomes constant after the acceleration returns to zero, contrary to the actual measured impulses (Yoganandan et al., 2008). In the case of kinematics-based injury metrics, the discrepancies in misidentifying concussive and subconcussive cases by HIC and RIC could be explained by these metrics' dependence on the shape of the acceleration impulse. In contrast, BrIC is only a function of peak angular velocity, and therefore exhibits less sensitivity to the shape of the impulse and the biphasic impulse's lack of restitutive deceleration (Figure 7).

In the case of brain angle metric, the biphasic impulse approximations show over five-fold errors compared to PCA-based impulses. This difference could be attributed to the simplification of the biphasic models that influences the solution of the mechanical lumped-parameter models. This discrepancy seems to affect the coronal direction the most and the sagittal direction the least for the biphasic approximations.

Similar to the brain angle results, brain finite element strains showed superior performance by our PCA-based approach. In previous publications, it has been shown that the closer the model is to the correct anatomical and morphological attributes of the brain, the more accurate the injury predictions become (Giordano et al., 2014; Zhao and Ji, 2019). Similarly, we see a decline in performance for the biphasic models as we include more region-specific and morphological details: from maximum principal strains (MPS) in the whole brain to MPS in corpus callosum and on to axonal fiber strains (Table 3).

Several other points are worth noting based on our analysis. Relative with time domain, we observe differences in kinematics profiles in the three anatomical directions. It seems that the head experiences higher linear accelerations in the anterior-posterior direction (3.54 ± 2.97 m/s²) compared to lateral (2.62 ± 2.25 m/s²) and inferior-superior (3.01 ± 3.44 m/s²) directions, and higher angular accelerations in the sagittal direction (190.53 ± 177.34 rad/s²) compared to coronal (314.50 ± 401.53 rad/s²) and axial (177.34 ± 128.16 rad/s²) directions. This observation might be attributed to the type and direction of loading in the specific activity, e.g., direction of tackling in football, as well as anatomical features such as the neck constraint in those directions (Eckersley et al., 2017). In the frequency domain, there is a dominant low-frequency response (10 to 40 Hz) in the head kinematics, expressed by ~ 90% of the total angular velocity response in the axial plane, ~ 83% in the coronal plane and ~74% in the sagittal plane, confirming previous findings on frequency dependence of head impacts (Wu et al., 2016; Laksari et al., 2018). These results could prove useful for designing better helmets and other safety devices to avoid brain injury by targeting specific low-frequency range of motion.

In summary, our proposed PCA decomposition approach not only provides a deeper understanding of the head's response during impacts, but also provides a formal basis for reconstructing and augmenting head impact kinematics data. Our current emulator is built upon the 537 measured on-field head impacts described above and successfully generates new kinematic data, whose features such peak or duration are nearly indistinguishable (see Supplementary Figure 7). It is expected that with more on-field measurements, we would be able to improve the performance of the emulator even further, but our convergence analysis showed the available 537 cases to be sufficient (see Supplementary Figure 8). This type of approach might prove necessary given the increased need for larger training data sets in modern machine learning algorithms.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

PA: study design, data analysis, and writing of the paper. KL: study design, data collection and analysis, and writing of the paper. HB and NT: study design and contribution to writing. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Prof. David Camarillo of Stanford University for providing the head kinematics. PA thanks the Fulbright and CONICYT project Equal Opportunities Scholarship 56170002, and also Carissa Grijalva of the University of Arizona, for helping with the visualizations. KL thanks BIO5 Institute at the University of Arizona for partial support for this study under award number 1119000-KL2.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2020.555493/full#supplementary-material

Footnotes

1. ^https://uweb.engr.arizona.edu/~klaksari/Resources.html

References

Abderezaei, J., Zhao, W., Grijalva, C. L., Fabris, G., Ji, S., Laksari, K., et al. (2019). Nonlinear dynamical behavior of the deep white matter during head impact. Phys. Rev. Appl. 12:014058. doi: 10.1103/PhysRevApplied.12.014058

CrossRef Full Text | Google Scholar

Arrué, P. (2020). Data-Driven Head Impacts Emulator. Laksari's Lab. Available online at: https://uweb.engr.arizona.edu/~klaksari/Resources.html

Babaee, H. (2019). An observation-driven time-dependent basis for a reduced description of transient stochastic systems. Proc. R. Soc. A 475:20190506. doi: 10.1098/rspa.2019.0506

PubMed Abstract | CrossRef Full Text | Google Scholar

Babaee, H., Choi, M., Sapsis, T. P., and Karniadakis, G. E. (2017). A robust bi-orthogonal/dynamically-orthogonal method using the covariance pseudo-inverse with application to stochastic flow problems. J. Comput. Phys. 344, 303–319. doi: 10.1016/j.jcp.2017.04.057

CrossRef Full Text | Google Scholar

Babaee, H., and Sapsis, T. P. (2016). A minimization principle for the description of modes associated with finite-time instabilities. Proc. R. Soc. Lond. A 472:2186. doi: 10.1098/rspa.2015.0779

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, Y., Wu, S., Zhao, W., Li, Z., Wu, Z., and Ji, S. (2018). Concussion classification via deep learning using whole-brain white matter fiber strains. PloS ONE 13:e0197992. doi: 10.1371/journal.pone.0197992

PubMed Abstract | CrossRef Full Text | Google Scholar

Daniel, W. (1990). Applied Nonparametric Statistics. Duxbury Advanced Series in Statistics and Decision Sciences. Boston, MA: PWS-KENT.

Google Scholar

Defense and Veterans Brain Injury Center (2018). DoD Worldwide Numbers for TBI | DVBIC.

DeKosky, S. T., Blennow, K., Ikonomovic, M. D., and Gandy, S. (2013). Acute and chronic traumatic encephalopathies: pathogenesis and biomarkers. Nat. Rev. Neurol. 9, 192–200. doi: 10.1038/nrneurol.2013.36

PubMed Abstract | CrossRef Full Text | Google Scholar

Eckersley, C. P., Nightingale, R. W., Luck, J. F., and Bass, C. R. (2017). “Effect of neck musculature on head kinematic response following blunt impact,” in Conference proceedings International Research Council on the Biomechanics of Injury, IRCOBI (Antwerp: IRCOBI).

Google Scholar

Gabler, L. F., Crandall, J. R., and Panzer, M. B. (2018a). Development of a metric for predicting brain strain responses using head kinematics. Ann. Biomed. Eng. 46, 972–985. doi: 10.1007/s10439-018-2015-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Gabler, L. F., Joodaki, H., Crandall, J. R., and Panzer, M. B. (2018b). Development of a single-degree-of-freedom mechanical model for predicting strain-based brain injury responses. J. Biomech. Eng. 140:031002. doi: 10.1115/1.4038357

PubMed Abstract | CrossRef Full Text | Google Scholar

Giordano, C., Cloots, R. J., van Dommelen, J. A., and Kleiven, S. (2014). The influence of anisotropy on brain injury prediction. J. Biomech. 47, 1052–1059. doi: 10.1016/j.jbiomech.2013.12.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Giordano, C., and Kleiven, S. (2014). Connecting fractional anisotropy from medical images with mechanical anisotropy of a hyperviscoelastic fibre-reinforced constitutive model for brain tissue. J. R. Soc. Interface 11:20130914. doi: 10.1098/rsif.2013.0914

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernandez, F., Wu, L. C., Yip, M. C., Laksari, K., Hoffman, A. R., Lopez, J. R., et al. (2014). Six degree of freedom measurements of human mild traumatic brain injury. Ann. Biomed. Eng. 43, 1918–1934. doi: 10.1007/s10439-014-1212-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, S., and Zhao, W. (2014). A pre-computed model responses atlas for instantaneous brain strain estimation in contact sports. Ann. Biomed. Eng. 43, 1877–1895. doi: 10.1007/s10439-014-1193-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, S., Zhao, W., Ford, J. C., Beckwith, J. G., Bolander, R. P., Greenwald, R. M., et al. (2014). Group-wise evaluation and comparison of white matter fiber strain and maximum principal strain in sports-related concussion. J. Neurotrauma 14, 1–43. doi: 10.1089/neu.2013.3268

PubMed Abstract | CrossRef Full Text | Google Scholar

Kimpara, H., and Iwamoto, M. (2012). Mild traumatic brain injury predictors based on angular accelerations during impacts. Ann. Biomed. Eng. 40, 114–126. doi: 10.1007/s10439-011-0414-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleiven, S. (2007). Predictors for traumatic brain injuries evaluated through accident reconstructions. Stapp Car Crash J. 51, 81–114. doi: 10.4271/2007-22-0003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleiven, S. (2013). Why most traumatic brain injuries are not caused by linear acceleration but skull fractures are. Front. Bioeng. Biotechnol. 1:15. doi: 10.3389/fbioe.2013.00015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kornhauser, M. (1954). Prediction and evaluation of sensitivity to transient accelerations. J. Appl. Mech. 21, 371–380.

Kuo, C., Wu, L. C., Ye, P. P., Laksari, K., Camarillo, D. B., and Kuhl, E. (2017). Pilot findings of brain displacements and deformations during roller coaster rides. J. Neurotrauma. 34, 3198–3205. doi: 10.1089/neu.2016.4893

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurt, M., Laksari, K., Kuo, C., Grant, G. A., and Camarillo, D. B. (2017). Modeling and optimization of airbag helmets for preventing head injuries in bicycling. Ann. Biomed. Eng. 45, 1148–1160. doi: 10.1007/s10439-016-1732-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Laituri, T. R., Henry, S., Pline, K., Li, G., Frankstein, M., and Weerappuli, P. (2016). New risk curves for NHTSA's brain injury criterion (BrIC): derivations and assessments. Stapp Car Crash J. 60, 301–362. doi: 10.4271/2016-22-0012

PubMed Abstract | CrossRef Full Text | Google Scholar

Laksari, K., Fanton, M., Wu, L., Nguyen, T., Kurt, M., Giordano, C., et al. (2019). Multi-directional dynamic model for traumatic brain injury detection. J. Neurotrauma 37, 982–993. doi: 10.1089/neu.2018.6340

PubMed Abstract | CrossRef Full Text | Google Scholar

Laksari, K., Kurt, M., Babaee, H., Kleiven, S., and Camarillo, D. (2018). Mechanistic insights into human brain impact dynamics through modal analysis. Phys. Rev. Lett. 120:138101. doi: 10.1103/PhysRevLett.120.138101

PubMed Abstract | CrossRef Full Text | Google Scholar

Laksari, K., Wu, L. C., Kurt, M., Kuo, C., and Camarillo, D. C. (2015). Resonance of human brain under head acceleration. J. R. Soc. Interface 12:20150331. doi: 10.1098/rsif.2015.0331

PubMed Abstract | CrossRef Full Text | Google Scholar

Low, T., and Stalnaker, R. (1987). “A lumped parameter approach to simulate the rotational head motion,” in International Research Council on Biokinetics of Impacts (IRCOBI), 203–215.

Google Scholar

Manoogian, S., McNeely, D., Duma, S., Brolinson, G., and Greenwald, R. (2006). “Head acceleration is less than 10 percent of helmet acceleration in football impacts,” in Biomedical Sciences Instrumentation.

PubMed Abstract | Google Scholar

Marjoux, D., Baumgartner, D., Deck, C., and Willinger, R. (2008). Head injury prediction capability of the HIC, HIP, SIMon and ULP criteria. Acc. Anal. Prevent. 40, 1135–1148. doi: 10.1016/j.aap.2007.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, L. E., Kuo, C., Wu, L. C., Urban, J. E., Camarillo, D. B., and Stitzel, J. D. (2018). Validation of a custom instrumented retainer form factor for measuring linear and angular head impact kinematics. J. Biomech. Eng. 140, 1–6. doi: 10.1115/1.4039165

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, L. E., Pinkerton, E. K., Fabian, K. C., Wu, L. C., Espeland, M. A., Lamond, L. C., et al. (2019). Characterizing head impact exposure in youth female soccer with a custom-instrumented mouthpiece. Res. Sports Med. 1–17. doi: 10.1080/15438627.2019.1590833

PubMed Abstract | CrossRef Full Text | Google Scholar

National Operating Committee on Standards for Athletic (2012). Standard Performance Specification for Newly Manufactured Football Helmets NOCSAE doc (ND) 002-11m12.

Newman, J. A., and Shewchenko, N. (2000). “A proposed new biomechanical head injury assessment function - the maximum power index,” in SAE Technical Papers. doi: 10.4271/2000-01-SC16

PubMed Abstract | CrossRef Full Text | Google Scholar

NTSA (1972). National Highway Traffic Safety Administration.

Google Scholar

Patil, P., and Babaee, H. (2020). Real-time reduced-order modeling of stochastic partial differential equations via time-dependent subspaces. J. Comput. Phys. 415:109511. doi: 10.1016/j.jcp.2020.109511

CrossRef Full Text | Google Scholar

Patton, D. A., McIntosh, A. S., Kleiven, S., and Fréchéde, B. (2012). Injury data from unhelmeted football head impacts evaluated against critical strain tolerance curves. Proc. Instit. Mech. Eng. 226, 177–184. doi: 10.1177/1754337112438305

PubMed Abstract | CrossRef Full Text | Google Scholar

Sapsis, T., and Lermusiaux, P. (2009). Dynamically orthogonal field equations for continuous stochastic dynamical systems. Phys. D 238, 2347–2360. doi: 10.1016/j.physd.2009.09.017

CrossRef Full Text | Google Scholar

Selassie, A. W., Wilson, D. A., Pickelsimer, E. E., Voronca, D. C., Williams, N. R., and Edwards, J. C. (2013). Incidence of sport-related traumatic brain injury and risk factors of severity: a population-based epidemiologic study. Ann. Epidemiol. 23, 750–756. doi: 10.1016/j.annepidem.2013.07.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegkas, P., Sharp, D. J., and Ghajari, M. (2019). The traumatic brain injury mitigation effects of a new viscoelastic add-on liner. Sci. Rep. 9:3471. doi: 10.1038/s41598-019-39953-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sullivan, S., Eucker, S. A., Gabrieli, D., Bradfield, C., Coats, B., Maltese, M. R., et al. (2015). White matter tract-oriented deformation predicts traumatic axonal brain injury and reveals rotational direction-specific vulnerabilities. Biomech. Model. Mechanobiol. 14, 877–896. doi: 10.1007/s10237-014-0643-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Takhounts, E. G., Craig, M. J., Moorhouse, K., McFadden, J., and Hasija, V. (2013). Development of brain injury criteria (BrIC). Stapp Car Crash J. 57, 243–266.

PubMed Abstract | Google Scholar

Taylor, C. A., Bell, J. M., Breiding, M. J., and Xu, L. (2017). Traumatic brain injury-related emergency department visits, hospitalizations, and deaths - United States, 2007 and 2013. MMWR Surveill. Summ. 66, 1–16. doi: 10.15585/mmwr.ss6609a1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, L. C., Kuo, C., Loza, J., Kurt, M., Laksari, K., Yanez, L. Z., et al. (2018). Detection of American football head impacts using biomechanical features and support vector machine classification. Sci. Rep. 8, 1–14. doi: 10.1038/s41598-017-17864-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, L. C., Nangia, V., Bui, K., Hammoor, B., Kurt, M., Hernandez, F., et al. (2016). In vivo evaluation of wearable head impact sensors. Ann. Biomed. Eng. 44, 1234–1245. doi: 10.1007/s10439-015-1423-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, L. C., Zarnescu, L., Nangia, V., Cam, B., and Camarillo, D. B. (2014). A head impact detection system using SVM classification and proximity sensing in an instrumented mouthguard. IEEE Trans. Bio-med. Eng. 61, 2659–2668. doi: 10.1109/TBME.2014.2320153

CrossRef Full Text | Google Scholar

Wu, S., Zhao, W., Ghazi, K., and Ji, S. (2019a). Convolutional neural network for efficient estimation of regional brain strains. Sci. Rep. 9:17326. doi: 10.1038/s41598-019-53551-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, S., Zhao, W., Rowson, B., Rowson, S., and Ji, S. (2019b). A network-based response-feature matrix as a brain injury metric. Biomech. Model. Mechanobiol. 19, 927–942. doi: 10.1007/s10237-019-01261-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoganandan, N., Li, J., Zhang, J., Pintar, F. A., and Gennarelli, T. A. (2008). Influence of angular acceleration-deceleration pulse shapes on regional brain strains. J. Biomech. 41, 2253–2262. doi: 10.1016/j.jbiomech.2008.04.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, F., Wu, Y., Norton, I., Rigolo, L., Rathi, Y., Makris, N., et al. (2018). An anatomically curated fiber clustering white matter atlas for consistent white matter tract parcellation across the lifespan. NeuroImage 179, 429–447. doi: 10.1016/j.neuroimage.2018.06.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, W., Ford, J. C., Flashman, L. A., McAllister, T. W., and Ji, S. (2016). White matter injury susceptibility via fiber strain evaluation using whole-brain tractography. J. Neurotrauma 33, 1834–1847. doi: 10.1089/neu.2015.4239

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, W., and Ji, S. (2019). White matter anisotropy for impact simulation and response sampling in traumatic brain injury. J. Neurotrauma 36, 250–263. doi: 10.1089/neu.2018.5634

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, W., Kuo, C., Wu, L., Camarillo, D. B., and Ji, S. (2017). Performance evaluation of a pre-computed brain response atlas in dummy head impacts. Ann. Biomed. Eng. 45, 2437–2450. doi: 10.1007/s10439-017-1888-3

PubMed Abstract | CrossRef Full Text | Google Scholar