Survival Model for Foot and Leg High Rate Axial Impact Injury Data

Objectives: Understanding how lower extremity injuries from automotive intrusion and underbody blast (UBB) differ is of key importance when determining whether automotive injury criteria can be applied to blast rate scenarios. This article provides a review of existing injury risk analyses and outlines an approach to improve injury prediction for an expanded range of loading rates. This analysis will address issues with existing injury risk functions including inaccuracies due to inertial and potential viscous resistance at higher loading rates. Methods: This survival analysis attempts to minimize these errors by considering injury location statistics and a predictor variable selection process dependent upon failure mechanisms of bone. Distribution of foot/ankle/leg injuries induced by axial impact loading at rates characteristic of UBB as well as automotive intrusion was studied and calcaneus injuries were found to be the most common injury; thus, footplate force was chosen as the main predictor variable because of its proximity to injury location to prevent inaccuracies associated with inertial differences due to loading rate. A survival analysis was then performed with age, sex, dorsiflexion angle, and mass as covariates. This statistical analysis uses data from previous axial postmortem human surrogate (PMHS) component leg tests to provide perspectives on how proximal boundary conditions and loading rate affect injury probability in the foot/ankle/leg (n = 82). Results: Tibia force-at-fracture proved to be up to 20% inaccurate in previous analyses because of viscous resistance and inertial effects within the data set used, suggesting that previous injury criteria are accurate only for specific rates of loading and boundary conditions. The statistical model presented in this article predicts 50% probability of injury for a plantar force of 10.2 kN for a 50th percentile male with a neutral ankle position. Force rate was found to be an insignificant covariate because of the limited range of loading rate differences within the data set; however, compensation for inertial effects caused by measuring the force-at-fracture in a location closer to expected injury location improved the model's predictive capabilities for the entire data set. Conclusions: This study provides better injury prediction capabilities for both automotive and blast rates because of reduced sensitivity to inertial effects and tibia–fibula load sharing. Further, a framework is provided for future injury criteria generation for high rate loading scenarios. This analysis also suggests key improvements to be made to existing anthropomorphic test device (ATD) lower extremities to provide accurate injury prediction for high rate applications such as UBB.


Introduction
Automobile safety regulations as well as design features in automobiles heavily rely on injury risk functions (IRFs) and associated injury criteria. Though multiple injury criteria may exist for the same body region, it is important to understand Associate Editor Clay Gabler oversaw the review of this article. Address correspondence to Ann M. Bailey, Department of Mechanical and Aerospace Engineering, University of Virginia, Center for Applied Biomechanics, 4040 Lewis and Clark Drive, Charlottesville, VA 22911. E-mail: amb9um@virginia.edu Color versions of one or more of the figures in the article can be found online at www.tandfonline.com/gcpi. how the variables chosen to predict injury relate to actual injury mechanisms since the International Organization for Standardization (ISO) does not specify a methodology for determining variables to include in the analysis (ISO 2014). Equally important is the ability of the anthropomorphic test device (ATD) to use these injury criteria to accurately interpret the injury risk for humans subjected to the same input conditions. This research focuses on leg fractures caused by axial impact loading. Though a bulk of leg injury research has focused on automotive intrusion, military underbody blast (UBB) produces injury to the same foot and leg bones (Alvarez 2011;Funk 2000), but it is not known whether existing injury criteria can be applied to these higher rate events. Automotive Table 1. Summary of previous injury risk functions Study Boundary condition Injury definition Injury predictor for 50% injury probability Yoganandan et al. (1999) Potted proximal tibia, unbooted, ballasted AIS 2+ 0.348 * Age + 0.415 * Axial tibia force = 4.4 6.8 kN axial mid-tibia force Funk et al. (2002) Potted proximal tibia, unbooted Bony foot/ankle fracture 45-year-old 50% male: 8.3 kN mid-tibia force 65-year-old 50% male: 6.1 kN mid-tibia force 45-year-old 5% female: 5.0 kN mid-tibia force 65-year-old 5% female: 3.7 kN mid-tibia force Bass et al. (2004) a Potted mid-femur, booted AFIS-S 2+ 8.6 kN mid-tibia force McKay and Bir (2009) a Potted femur, unbooted AFIS-S 4+ 6.4 kN axial tibia force, 10.8 m/s velocity Quenneville et al. (2010) Potted tibia, unbooted No foot Bony tibia fracture 7.9 kN = 10% risk of injury Henderson et al. (2013) a Potted proximal tibia, unbooted, ballasted AIS 2+ 7.34 kN distal tibia force 6.16 kN proximal tibia force Gallenberger et al. (2013) Potted proximal tibia, unbooted, ballasted Bony fracture Neutral position: 6.8 kN proximal tibia force 20 • dorsiflexion: 7.9 kN prox. tibia force Yoganandan et al. (2014) Potted proximal tibia, unbooted, ballasted AIS 2+ 25-year-old: 10.4 kN proximal tibia force 45-year-old: 8.3 kN proximal tibia force 65-year-old: 6.6 kN proximal tibia force a Studies aimed at replicating UBB or landmine blast events.
intrusion events normally yield floor accelerations less than 250 g with 10-50 ms durations and floor pan intrusions in excess of 350 mm (Crandall et al. 1998;Funk et al. 2002), whereas UBB events exceed 1,000 g of floor acceleration with durations as small as 3 ms and intrusions as low as 30 mm (Scherer et al. 2010). Several studies have determined IRFs based on various parameters measured in postmortem human surrogate (PMHS) experiments Table 1. Though these statistical models fit their respective injury data, their boundary conditions and loading rates should be considered before applying them to different applications. The purpose of this work is to generate an IRF for the leg that may be utilized for a broader range of applications rather than limited to the boundary conditions under which they were created. Evaluating several previously developed IRFs from a variety of loading rates, it is not clear which are applicable for the higher loading rates associated with UBB. Though the loading rates from Henderson et al. (2013) and McKay and Bir (2009) were among the highest of those compared in Figure 1, the location of force measurement is proximal to the location of a majority of the injuries produced. Injury risk functions based on mid-tibia forces measured by implanted tibia load cells reveal differences in injury tolerance between the McKay IRF developed at high load rates and the Funk et al. (2002) data from more moderate load rates. It is unclear whether these differences are due to input pulse characteristics, quality of specimens tested, or the way injury is defined in the 2 models. It is important to understand the effects of these dissimilarities in order to make progress in understanding the role of loading rate in leg injury.
One of the most widely used IRFs for the leg was generated by Yoganandan et al. (1996). Despite a relatively large sample size (n = 52), the data set is composed of force measurements taken from various locations and includes both dorsiflexed and neutral ankle positions. Dorsiflexion has since been shown to increase the fracture tolerance of the ankle joint (Klopp et al. 1997). Furthermore, the proximal boundary conditions for the included data were different, which enhances the potential for inertial differences between the tests and increases the need for a common location of force measurement. Though the data from Yoganandan et al. (1996) were reanalyzed to use distal forces (Yoganandan et al. 2014), the analysis still fails to account for the effect of ankle flexion.
Boundary conditions in tests aimed at determining injury risk are of vital importance; however, simulating exact boundary conditions and achieving the higher loading rates required for studying UBB is challenging. McKay andBir (2009), Funk et al. (2002), and Henderson et al. (2013) each achieved higher loading rates with 3 different boundary conditions, thus producing different injury patterns in select cases. All but Henderson et al. produced proximal tibia fractures, which was likely due to the presence of the knee joint in the other test series. Knowing the effect that each of these complicated boundaries has on producing unrealistic stresses on potentially weaker bony regions is important for determining whether these moreproximal injuries should be considered for inclusion in an analysis focusing on injuries to the foot/ankle complex.  Yoganandan et al. (1996) 16 kg ballast, potted proximal tibia 1.8-6.9 0.05-8.1 1.0-11.0 15 2 Klopp et al. (1997) Seated position a 1.4-5.5 0.2-7.2 1.6-10.8 15 10 Funk et al. (2000) Potted mid-femur a 4.8-7.1 1.1-4.7 4.4-10.0 2 15 Funk et al. (2002) Knee joint against fixed boundary 1.4-8.4 0.1-6.3 1.7-10.8 7 7 a Includes dorsi/plantar flexion tests.
McKay and Bir (2009) proposed a velocity-based injury metric for the leg to be used for UBB loading conditions. Though impact velocity is significantly correlated with injury in their data, the velocity-based injury criterion should not be used when boundary conditions vary or with different impactor masses. Velocity, as an injury predictor, neglects the effect of proximal boundary conditions on loading rate, and no validation has been performed to determine whether the proximal boundary conditions used in that study are realistic for a human in a UBB event. Using engineering principles to choose predictor variables can enhance the versatility of the IRF by enabling it to be used for loading scenarios other than those similar to the experimental configuration.
Some IRFs have taken into account covariates such as age, sex, and bone mineral density (Funk et al. 2002;Kitagawa et al. 1998;Klopp et al. 1997;Roberts 1993;Yoganandan et al. 1996). Though these works have used realistic injury predictors (force) and have attempted to compensate for specimen variability, the location of force measurement limits their accuracy as well as the relationship established between the main injury predictor and the covariates. For example, Funk et al. measure the force at mid-tibia using an implanted load cell. Though this force is realistic for determining force-at-fracture for the tibia and locations in close proximity to the mid-tibia, most injuries occurred at the level of the calcaneus. Further analysis of the Funk et al. data shows that footplate force is better correlated with injury than mid-tibia force. Load sharing between the tibia and fibula can account for about 6% of the difference between the load measured at the two locations (Funk et al. 2004;Goh et al. 1992;Takebe et al. 1984), and the remainder can likely be assigned to the effect of inertia due to the mass present between the 2 load cells. Table 2 show that the studies aimed at understanding the effect of UBB on the lower leg have different proximal boundary conditions and instrumentation locations than others performed at automotive rates, making it difficult to understand the effect of higher loading rates on the lower leg. As the rate of loading increases, it becomes increasingly important to measure the load closer to the region of interest because of inertial effects. As the acceleration of the leg increases, the greater the difference between the forces measured in one location versus another due to the mass present between the load cell and region of interest. Though previous IRFs may represent the data set from which they are generated, the predictive ability of that function for another data set decreases with deviations from the loading rate or boundary conditions. Thus, it is important to consider the methods and conditions under which an IRF was developed when applying it to a new data set and attempting to predict injury, particularly for the case of applying automotive rate injury risk curves to higher rate UBB conditions. The aim of this article is to provide a leg injury criterion using a physically reasonable predictor and covariates that can be applied to measurements from a biofidelic ATD exposed to various boundary conditions and load rates, while presenting a methodology that can be employed for other biomechanics applications. This analysis provides an injury criteria based on plantar force rather than mid-or proximal tibia force to decrease errors due to inertial effects, thus enhancing its accuracy for higher loading rates such as UBB.

Methods
Data were collected from 150 previous axial loading component PMHS leg tests (Funk et al. 2002;Gallenberger et al. 2013;Henderson et al. 2013;Klopp et al. 1997;Smith et al. 2005;Yoganandan et al. 1996), and injuries were evaluated to determine where fracture force should be measured. Injuries included mainly calcaneus, talus, and distal tibia (pilon) fractures, with calcaneus fractures far outnumbering other injuries. Further analysis showed that calcaneus fractures often occurred alone (38 tests) and also accompanied talus and pilon fractures. Proximal tibia injuries occurred only in tests in which the knee joint was included (Funk et al. 2002;McKay and Bir 2009) and occurred only once without a calcaneus fracture (Figure 2). Thus, plantar force, measured by footplate load cells, was chosen as the most sensible factor for predicting injury, due to the majority of bony fractures being located in the calcaneus. Further justification for choosing plantar force over mid-or proximal tibia force is provided by Funk et al. (2002), who stated that "[acoustic emission] data showed that calcaneal fracture correlated more with peak axial footplate force" (p. 756).
The data set was then reduced by eliminating tests that did not report footplate or plantar force or for which the injury mechanism was not axial loading. For example, Roberts (1993) reported injuries caused by ankle rotation due to excessive intrusion, which were removed from the data set. Dorsiflexion and plantar flexion tests were included to increase the sample size, and though injury tolerance has been shown to change due to differences in subtalar joint contact area, injuries from these positions cause injury to the same bony structures of the  Klopp et al. (1997), Smith et al. (2005), and Yoganandan et al. (1996); 150 tests.
foot as axial loading to a neutrally positioned ankle (Calhoun et al. 1994;Wang et al. 1995). The remaining 82 tests that were used in the survival analysis are summarized in Table 2 and  Table A1 (see online supplement).
In most studies, the Abbreviated Injury Scale (AIS) is used to define the division between "injury" and "no injury" (Gennarelli and Wodzin 2005;ISO 2014). Some applications for IRFs, however, may dictate division between acute and prolonged injury effects; for example, the distinction between combat injuries that prevent warfighters from escaping further threats after an initial injury event as in the case of UBB or anti-personnel landmines. Consideration must also be given to injury mechanism to ensure that injuries are distinguishable using a predictor variable or covariate. For the foot/ankle/leg complex, it is not reasonable to generate separate injury criteria for specific bones because of overlap in the injury thresholds for each bone (Kuppa et al. 2001). For example, for a given axial force, a talus injury might occur in one specimen while a calcaneus fracture might occur in another due to interspecimen variability. Further, development of a survival model should depend on a causal relationship between predictor variables and injury outcome or injury level (McMurry and Poplin 2014); thus, any bony fracture in the foot/ankle/leg was considered as "injury" for this analysis. The Ankle and Foot Injury Scale for Severity (AFIS-S) (Levine et al. 1995;Manoli et al. 1997) severity scores are shown in Table A1 to provide further information about the level of injury for each case.
Covariates were chosen with the goal of reducing variance in the relationship between injury and plantar force (Hougaard 2000). Before adding covariates to the model, it was important to understand how factors relate to each other as well as how they physically relate to the predictor variable and injury. Simply assuming that each covariate has a linear effect on injury outcome could greatly affect the predictive capability of the injury criterion. For this analysis, age (years), sex (male = 1, female = 0), whole-body mass (kg), a sex-mass interaction term, and dorsiflexion angle (degrees), an age-sex interaction term, and force rate were considered as variables affecting force-at-fracture. Flexion angle was included due to previous studies demonstrating its ability to increase force-atfracture (Gallenberger et al. 2013;Klopp et al. 1997); dorsiflexion was treated as positive, whereas plantar flexion angles were assigned negative values. Force rate was not included because of its small contribution to the yield strength of bone under the range of loading rates included in the data set (McElhaney 1966). Use of plantar foot force rather than proximal or mid-tibia force allowed for a less rate-dependent force measurement because of smaller contributions from inertial differences between load measurement location and injury location. Further, mass to the two-thirds power was used as a covariate to account for differences in force measured in different sized specimens because of mass's dimensional relationship to force (Eppinger et al. 1984). Mass-normalized plantar force was considered in place of mass as a covariate; however, the ability to account for mass separately in an injury criterion was considered advantageous to the model's versatility, particularly when determining the effect on separate sexes. Though the Yoganandan et al. (1996) tests attempted to mass-normalize by ballasting the specimens to 16 kg, the mass covariate in this case is meant to compensate for individual leg characteristics that impact how specimen mass and mass recruitment affect fracture force. The interaction term for sex and mass was considered for this reason. Age, sex, mass to the two-thirds power, and dorsiflexion angle were included as covariates to plantar force in the selected survival model.
Statistical software R version 3.1.1 (R Foundation, Vienna, Austria) was used to perform the survival analysis using a multivariate Weibull distribution. The Akaike information criterion (AIC) was used to compare potential survival distributions in accordance with recommendations by the ISO (2014). A Weibull distribution (AIC = 473.53) was chosen over lognormal (AIC = 475.80) and log-logistic (AIC = 476.89) distributions; however, it should be noted McMurry and Poplin (2014) demonstrated through simulation that AIC makes a somewhat arbitrary choice with typical sample sizes. Nonetheless, the selected Weibull model is consistent with their recommendations associated with "aesthetic considerations" (p. 624). Plantar force was used as the main predictor variable. Tests in which the specimen was not injured were considered right-censored (39 cases), whereas tests in which the specimen was injured but the exact force-at-fracture is unknown were treated as left-censored (22 cases). The 21 tests for which force-at-fracture was known were treated as uncensored.

Results
The studies selected for the survival analysis consisted of axial impact loading to the leg with neutral or flexed ankle positioning and minimal rotation. The data set consisted of 37 females with an average age (±standard deviation) of 70 ± 11 years and average body mass of 59.5 ± 5.1 kg, whereas the 45 males had an average age of 60 ± 14 years and average mass of 78.2 ± 13.1 kg. A summary of the correlations between predictor variables considered in this analysis is presented in Table A2 (see online supplement), and Figure A1 (see online supplement) shows the distribution of age, dorsiflexion angle, footplate forces, and rates within the data set. Sex and mass were correlated within the data set (r = 0.675) and age and sex were slightly correlated (r = −0.349).
Parameters for the selected model are provided in Table 3 along with P values and standard errors for each model coefficient to be used with Eqs. (1)-(3). In Eq. (1), F represents plantar force (the main predictor variable), and β i represent the coefficients estimated for the survival model, and P i in Eq.
(2) represent covariates chosen for the analysis. Dorsiflexion was found to be significant to a level of α = .05, and age was considered to be marginally significant (P = .072). Sex and mass were not deemed statistically significant but remain in the model to adjust for any unmeasured potential confounding effects; these variables have a moderate effect size based on correlations with plantar force shown in Table A2, which indicates potential for significance provided additional data to increase the power. Each model covariate behaved as expected within the model. Male specimens of greater mass and larger dorsiflexion angles are associated with an increase the force-at-fracture for a given probability of injury, whereas age decreases the fracture force. The age-sex interaction term was found to be insignificant (P = .63) and therefore was not included in the final model. lso presents tabulated forces for 10, 50, and 90% risk of injury for different age, sex, mass, and dorsiflexion angle values. Figure 3 shows a comparison of the IRF for 45-year-old males and females, as well as a comparison of how dorsiflexion angle affects the probability of injury.

Discussion
Comparison to similar IRFs reveals a higher predicted fracture force for a 50% probability of injury, as expected because of proximity to the location of impact. For a 50th percentile male, the current statistical model predicts a 50% probability of injury given a plantar force of 10.2 kN, compared to 8.3 kN predicted by Funk et al. (2002) and Yoganandan et al. (2014) at the mid-and proximal tibia ( Figure 4). These results are reasonable when considering that inertia due to location of the load cell could add, on average, 2.6 kN to the mid-tibia force, assuming the mass inferior to the mid-tibia is 1.96 kg for a 75 kg male (one quarter of the mass of the leg plus the foot mass; Plagenhoef et al. 1983), and the average acceleration at the mid-tibia is around 136 g (calculated from Funk et al.'s [2002] accelerations). The mid-tibia force predicted by McKay and Bir (2009) for 50% probability was even lower (6.4 kN); though tibia accelerations were not provided for that study, the presence of proximal tibia fractures as well as potentially higher accelerations in the tibia could lead to an even greater difference between plantar and mid-tibia forces. This justification may also explain the lack of statistical significance of force as an injury predictor in the McKay and Bir data set, which spans velocities of 7-12 m/s. Differences in 50% probability of injury force levels at the mid-tibia further warrant an injury predictor that is less influenced by loading rate and  boundary conditions. The survival model presented here overcomes these obstacles, which are some of the major divisions between UBB and automotive intrusion applications. An age-sex interaction term was not used in this analysis due to the model predicting an increase in fracture force for older women, which is not consistent with the expected outcome and previous studies (Keaveny et al. 2009;Launey et al. 2010). Instead, by excluding the age-sex interaction term, the model assumes that fracture tolerance varies the same with age for men and women. This assumption is likely reasonable considering that the sample population has excluded osteoporotic specimens, predisposing the sample to exclude the postmenopausal females, who cause a divergence in the bone strength-age relationship between men and women.
A possible source of error within this survival analysis is the use of whole-body mass rather than component leg mass as a covariate. Though whole-body mass allows for some correction for size/mass within the analysis, it does not account for distribution of the mass within the body, which could lead to over-or underestimation of specimen mass when trying to scale the force results. Additionally, mass distribution is different between men and women; thus, with more available data, using a sex-mass interaction term would be appropriate, particularly if mass and sex within the population are correlated. Further fracture force inaccuracy may arise from injury occurring in more proximal to the calcaneus. Future work will aim at determining the magnitude of this error.
The major contribution of this study is a statistical model that can be used to estimate probability of injury to the foot/ankle/leg due to axial loading with decreased sensitivity to loading rate or boundary conditions, which are the main contributors to differences between UBB and automotive intrusion scenarios. Though additional uncensored data from a more representative sample of the population would improve the predictive properties of the model, such data are expensive and difficult to collect using PMHS testing. Implications of the presented survival model include a need for an improved ATD leg that incorporates accurate transmissibility of force over a range of loading rates, as well as a load cell located on the plantar surface of the foot. Though previous studies were based on mid-tibia force primarily for applying the injury criteria to ATDs, modification of the current ATD legs is recommended to allow for use of this more comprehensive IRF. The proposed IRF would allow for more accurate assessment of injury mitigation using ATDs, specifically because of its validity at both automotive and military UBB rates of loading. Forthcoming research efforts will include characterizing how the presented IRF behaves when used with existing ATD legs, in addition to an evaluation of the model's predictive capabilities by means of future PMHS injury data.