Integrating Test Data and ATBM Simulations into Dose Propagation Uncertainty Formulation for Bone Fracture Risk Assessment

Hongyun Wang¹, Corinne Kramer², Jessica Swallow², Wesley A. Burgei³, Hong Zhou^4*

¹Department of Applied Mathematics, University of California, Santa Cruz, CA, USA.
²Institute for Defense Analyses, Alexandria, VA, USA.
³U.S. Department of Defense, Joint Non-Lethal Weapons Directorate, Quantico, VA, USA.
⁴Department of Applied Mathematics, Naval Postgraduate School, Monterey, CA, USA.

DOI: 10.4236/health.2019.1110109                

Abstract

We consider the problem of assessing bone fracture risk for a subject hit by a blunt impact projectile. We aim at constructing a framework for integrating test data and Advanced Total Body Model (ATBM) simulations into the risk assessment. The ATBM is a finite element model managed by the Joint Non-Lethal Weapons Directorate for the purpose of assessing the risk of injury caused by blunt impacts from non-lethal weapons. In ATBM simulations, the quantity that determines arm bone fracture is the calculated maximum strain in the bone. The main obstacle to accurate prediction is that the calculated strain is incompatible with the measured strain. The fracture strain measured in bending tests of real bones is affected by random inhomogeneity in bones and uncertainty in measurement gauge attachment location/orientation. In contrast, the strain calculated in ATBM simulations is based on the assumption that all bones are perfectly elastic with homogeneous material properties and no measurement uncertainty. To connect test data and ATBM simulations in a proper and meaningful setting, we introduce the concept of elasticity-homogenized strain. We interpret test data in terms of the homogenized strain, and build an empirical dose-injury model with the homogenized strain as the input dose for predicting injury. The maximum strain calculated by ATBM has randomness due to uncertainty in specifications of ATBM setup parameters. The dose propagation uncertainty formulation accommodates this uncertainty efficiently by simply updating the shape parameters in the dose-injury model, avoiding the high computational cost of sampling this uncertainty via multiple ATBM runs.

Keywords

Blunt Impact, Bone Fracture, Effects of Loading Rate, Elasticity-Homogenized Strain, Dose Propagation Uncertainty Formulation, ATBM Simulations

1. Introduction

Non-lethal blunt impact weapons have been widely used by law enforcement and military to incapacitate individuals while minimizing fatalities and collateral damage.

With continued use of existing blunt impact projectiles and the development of new capabilities, it is important to assess risk of injury. The Joint Non-Lethal Weapons Directorate (JNLWD) has been developing blunt impact injury modeling capabilities for over a decade beginning with the Interim Total Body Model (ITBM). ITBM eventually evolved into the Advanced Total Body Model (ATBM), which is a simulation package still in development today. The existing simulation engine consists of two parts: a finite element model (FEM) based on a 50th percentile male target; with individual tissue models and properties defined from the literature. This model is used in conjunction with a projectile model to simulate the impact dynamics including forces, deformation, stresses and strains in individual elements over time. The second part is called the injury model, usually a logistic regression of experimental injury data [1]. A probability of hit model is also under development for this package. The ATBM package includes models for the head, eye, neck, thorax, arm and leg. The output values of a given simulation are compared to the injury model to assess the probability of injury under the simulated conditions. We note two limitations to this methodology, which together, result in significant uncertainty in the predicted risk of injury. They arise from the fact that the experimental population is both limited in quantity (small N) and different from the desired target population.

The FEM consists of body part properties based on the best available data and 50th percentile geometry. All estimates of response metrics are intended to represent the 50th percentile response; however we have limited data on how such properties are distributed across the population, and we do not have adequate confidence in the value of the mean (or median) property. This is a fundamental limitation that can only be mitigated through increased experimental efforts, which come at significant cost. In addition, the properties that are available are mainly sourced from post-mortem human subject (PMHS) experiments, which should be assumed to have some differences from living subjects. So the FEM must extrapolate from already highly uncertain values.

The injury model is developed using all available data, some of which is in the correct strain rate regime¹ and all of which is obtained on PMHS or animal subjects. Again, by necessity, we extrapolate our injury model results from the experimental population to the target.

In this study, we look at the issue of uncertainty and extrapolation in the injury model, as a venue for connecting the idealized ATBM simulations and real test data. Based on the literature, we expect that skeletal injuries are likely better characterized by the experimental data set than are soft tissue injuries (due to the fact that bone properties are better maintained through death and preservation than soft tissues). Therefore here we focus on predicting bone fracture risk, and we examine uncertainty propagation. Recently, Kramer and Swallow have raised the issue of sensitivity in the ATBM simulation outputs with respect to the inputs. They explored risk assessment of significant injury using ATBM simulations and studied the propagation from uncertainty in ATBM simulation results to the predicted distribution of injury risk and the average risk (unpublished results). In this study, we aim at establishing a mathematical framework for assessing bone fracture risk that incorporates theoretical ATBM simulations and experimental test data of real bones. Specifically, we study arm fracture caused by the impact of non-lethal weapon projectiles. We build on previous work [2] that endeavors to incorporate uncertainty into the estimates of fracture risk. We also examine several mathematical issues in combining ATBM simulations and test data.

The experimental data on arm fracture were measured in three-point bending tests on a small number of humeri and forearms from female cadaveric donors of average age 57 (humerus samples) and 61 (forearm samples) [3] [4]. The PMHS samples were tested under loading conditions that were designed to emulate an airbag impact in traffic accidents [3] [4] [5].

We model the data using a log-normal distribution of stresses, strains and bending moments and apply a linear correction for the documented gross physical characteristics of the subjects (age and body mass). This gives us a venue to correct for the effects of known characteristics in test data and to explore the variations in macroscopic properties that are not fully correlated with the age and body mass.

We explore the issues of including and adjusting for the effects of loading rate, age and body mass in the risk assessment framework. The goal is to formulate an injury risk prediction that is relevant and meaningful for the larger portion of the population. The key step in integrating test data and ATBM simulations is to recognize that the calculated strain and the measured strain are different. The fracture strain measured in drop tests of real bones is affected by random inhomogeneity in bone and uncertainty in measurement gauge attachment (location/orientation). The strain calculated in ATBM simulations is for an idealized situation where all bones are perfectly elastic with homogeneous material properties and no measurement uncertainty. To assimilate test data with ATBM simulations, we introduce the concept of elasticity-homogenized strain, which is the hypothetical strain in experiments if the bone remains perfectly elastic with homogeneous material properties all the way up to the fracture. The elasticity-homogenized strain serves as a bridge connecting test data to ATBM simulations. The connection process consists of several steps: 1) interpret test data in terms of elasticity-homogenized strain; 2) in the framework of dose propagation uncertainty, build a dose-injury model with ln(elasticity-homogenized strain) as the input dose for predicting injury; 3) in ATBM simulations, calculate the input dose as ln(maximum strain); and 4) predict the injury probability using the properly built dose-injury model and the output of ATBM simulations. The result is a computational framework for assessing the injury risk of bone fracture for a subject hit by a non-lethal weapon projectile.

We organize the paper as follows. First, we briefly review the setup of threepoint bending tests in Section 2, which provides the experimental data on impact-induced bone fracture of real arm samples with bones surrounded by tissues and muscles. In Section 3, we describe the general dose-response relation based on the uncertainty propagation formulation. The static dose response relation is applicable in assessing injury risk of bone fracture when the impact loading is quasi-static and the injury process is nearly memoryless. In a dynamic loading test, the test result is in the form of measured peak value of dose (strain or stress or bending moment) at fracture, which is called fracture tolerance. In Section 4, we develop an efficient and accurate inference procedure for determining model parameters in the dose-response relation from observed samples of fracture tolerance in loading tests. When the impact occurs at a sufficiently high loading rate, the injury process is no longer memoryless due to the dynamic plastic deformation of bones. Consequently, the fracture tolerance is affected by the loading rate. In Section 5, we derive scaling laws of fracture tolerance versus loading rate and we examine the scaling laws using data of compression tests conducted on human femoral cortical bone over a wide range of loading rates. In Section 6 and Section 7, we carry out fluctuation analysis on the existing data of humerus fracture and forearm fracture, and build empirical injury models that are compatible with the ATBM simulation environment. Section 8 synthesizes the empirical injury models and the ATBM simulations to produce a viable framework for assessing the bone fracture risk caused by the impact of a blunt projectile. Finally, we highlight our results and outline future work in Section 9.

2. Review of Three-Point Bending Test Setup

In three-point bending experiments, a rod-shaped sample is placed horizontally, supported/anchored at two ends. The load is applied in the vertical direction, at a location between the two ends. In [5] [6], dynamic three-point bending tests were performed on PMHS forearms by dropping a heavy impactor (9.48 kg) onto the anchored sample, from a specified height (2 m) above the sample. The drop height was selected to produce a moderate pre-impact velocity (3.63 or 4.42 m/s) for the impactor. If the sample is perfectly elastic for the full range of stress, upon the impact, it will deform elastically and when the load is released it will relax back to its original state with no damage. The forearm is not perfectly elastic. When the load is large enough, it will fracture. In the three-point bending tests of forearms [5] [6], upon the impact, the strain, stress and bending moment increase initially before the occurrence of fracture; they reach their peak values at fracture, and then decrease after fracture since the fractured sample can no longer sustain the load. In experiments, the mass and the drop height of impactor are selected to ensure that fracture will occur in the process; the peak values of strain, stress, and bending moment observed at fracture are recorded as fracture tolerances. The data set from three-point bending experiments is a sequence of random values of fracture tolerance, each corresponding to an individual bone sample tested. In injury risk assessment, a relevant task is to estimate the probability of bone fracture at a given input dose in the form of strain, stress, or bending moment. This task is different from the question directly addressed in the dynamic loading tests. Below we first introduce the mathematical formulation of dose injury relation, which predicts injury probability from the given input dose. Then we use bending test data to determine model parameters in the dose injury relation.

3. The Dose Injury Relation

3.1. Dose Quantity for Predicting Injury Risk

There are several candidates that may serve as the single metric predictor variable for injury risk assessment. In the three-point bending tests, these candidates include strain, stress and bending moment. Any of these predictor variables or any function of these variables may be considered as the input dose in the mathematical formulation for predicting the probability of bone fracture injury. In particular, we can use ln(strain) instead of strain as the input dose. Mathematically, let

x = quantity selected as the input dose.

Here the phrase “input” refers to the situation where the dose x is specified and maintained in experiments. For example, in a clinic trial of a drug, the input dose is the daily amount of medicine administrated to each patient, which is precisely specified and maintained. In some experiments, the dose x is not directly specified. Instead, a related quantity u is directly specified and quantity x is a consequence of quantity u and is measured. Consider a hypothetical test setup for investigating the relationship between the bending moment and the fracture probability. The impactor is dropped from various vertical positions, spanning a wide range of height. In each repeat of the hypothetical test, the drop height is prescribed, and the corresponding maximum bending moment on the sample is measured regardless of the injury outcome. The maximum bending moment during each loading serves as the dose for predicting fracture injury. The bending moment is positively correlated with the prescribed drop height of impactor, but is not completely determined by the height. When the dose quantity x is not directly specified but is influenced by a specified quantity, it is more appropriate to label quantity x as “measured” dose. For simplicity of discussion, most of the time, we shall use these two terms interchangeably. The data set of the hypothetical test has the form

$\text{Data}=\left\{x_j,I_j,j=1,2,\cdots ,m\right\}$

where

m = number of trials repeated in experiments;

x_j = measured (or specified) dose in the j-th trial;

I_j = observed binary injury outcome in the j-th trial.

In this hypothetical test setup for measuring the injury probability, to capture the variation of injury probability vs dose, it is desirable to select a wide range of drop heights to cover the section from near 0% to near 100% injury probability.

Notice that the three-point bending test in [5] [6] is designed quite differently from the hypothetical test above. In [5] [6], the drop height of the impactor is set high enough to ensure fracture occurs in each and every loading test. The dose increases monotonically until fracture, and the peak dose at fracture is recorded as the fracture tolerance. Such an experimental setup is called dynamic loading. A key characteristic of dynamic loading is the virtually unlimited increase of load to ensure eventual fracture in each test. Accordingly, the data set measured in dynamic loading experiments takes the form

$\text{Data}=\left\{Y_j,j=1,2,\cdots ,m\right\}$

where

m = number of samples tested;

Y_j = observed fracture tolerance of sample j;

= measured value of dose x at fracture of sample j.

Here we use x to denote the specified/measured dose and Y to denote the observed fracture tolerance (i.e., the value of x at fracture). This distinction in notation is necessary for interpreting the results of dynamic bending tests.

3.2. Logistic Injury Model

Let p(x) be the probability governing the occurrence of injury at dose quantity x.

We model the injury probability using the logistic function [2] :

$p\left(x\right)=\frac{1}{1+\exp \left(\frac{-2\ln \left(9\right)}{W}\left(x-D_{50}\right)\right)}\equiv f_L\left(x;D_{50},W\right)$ (1)

where

D₅₀ = the median injury dose, at which injury probability p(D₅₀) = 50%;

in general, D_η = the dose value at which p(D_η) = η%; and

W ≡ D₉₀ – D₁₀ = the 10 - 90 percentile width of the injury function.

3.3. Formulation of Dose Propagation Uncertainty

In a previous study [2], we interpreted the logistic injury model as a very accurate approximation to the normal distribution model, which is formulated based on: 1) uncertainty in dose propagation; 2) uncertainty in dose measurement; and 3) uncertainty in critical dose threshold for injury. The formulation of dose propagation uncertainty is summarized as follows.

· The target dose is the amount of dose that propagates from the input site to reach the active site for injury. At the active site, the target doze Z and the critical threshold Z^(c) uniquely determines the binary injury outcome I:

$I=\{\begin{array}{l}1\left(\text{injured}\right),\text{if}Z>Z^{\left(c\right)}\\ 0\left(\text{notinjured}\right),\text{otherwise}\end{array}$

Once the target dose Z and the critical threshold Z^(c) are known, the injury outcome is completely determined.

· The critical threshold Z^(c) is not a fixed value for all samples. Even for bone samples from donors virtually of the same age, of the same bone density and of the same bone geometry, random microstructural constituents of samples will result in different values for the critical threshold. At a fixed target dose Z, even for a set of seemingly homogeneous samples, the hidden variation in the critical threshold Z^(c) will lead to different injury outcomes. For samples from a heterogeneous population, there will be additional variations in the critical threshold Z^(c), attributed to differences in bone density and bone geometry. We model Z^(c) as a Gaussian random variable.

$Z^{\left(c\right)}={\mu }_1+s_1{\xi }_1,{\xi }_1~N\left(0,1\right)$

· There is a discrepancy between the measured dose x at the measurement site and the target dose Z reaching the active site of injury. The discrepancy (Z – x) reflects:

1) The uncertainty in dose propagation from the measurement site to the active site for injury; and

2) The uncertainty in measurement gauge attachment (position and orientation).

We model the discrepancy (Z – x) as an independent Gaussian random variable

$\left(Z-x\right)=-{\mu }_2+s_2{\xi }_2,{\xi }_2~N\left(0,1\right)$

where (ξ₁, ξ₂) are i.i.d. samples of N(0, 1).

· The algebraic difference (Z – Z^(c)) has the expression

$Z-Z^{\left(c\right)}=x-\left({\mu }_1+{\mu }_2\right)-\left(s_1{\xi }_1-s_2{\xi }_2\right)\equiv \left(x-\mu \right)-s\xi$

where

$\mu =\left({\mu }_1+{\mu }_2\right),s=\sqrt{s_1^2+s_2^2},\xi ~N\left(0,1\right)$

That is, the difference (Z – Z^(c)) is another Gaussian random variable.

· When the measured dose is at level x, the injury probability of a random sample in a random realization of test, p(x), has the expression:

$p\left(x\right)=\mathrm{Pr}\left(Z-Z^{\left(c\right)}>0\right)=\mathrm{Pr}\left(\xi <\frac{x-\mu }{s}\right)=\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{x-\mu }{\sqrt{2}s}\right)$

· The injury model is completely specified by two shape parameters: median injury dose D₅₀ and the width W ≡ D₉₀ – D₁₀. It has the expression:

$p\left(x\right)=\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{2{\text{erf}}^{-1}\left(\text{0}\text{.8}\right)\cdot \left(x-D_{50}\right)}{W}\right)\equiv f_N\left(x;D_{50},W\right)$ (2)

where the shape parameters (D₅₀, W) are related to (μ, s) as

$D_{50}=\mu ,W=2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right)s$ (3)

Both the logistic model and the normal distribution model are completely specified by shape parameters (D₅₀, W). In [2], we showed that these two models are very close to each other. In the next section, we examine two inference approaches based on these two models.

4. Inferring the Dose-Injury Relation in a Memoryless Process

4.1. Memoryless Process

We consider the case where the applied dose varies with time but the dose injury process is still memoryless. Specifically, by “memoryless process”, we assume:

· The dose at time t takes effect immediately without delay;

· If the dose at time t does not result in injury immediately, it will not influence the injury process beyond time t; that is, there is no partial injury or left-over effect;

· The randomness in the occurrence of injury is associated with biovariability of the test sample and realization of the experimental setup, both of which are invariant with respect to time during each loading process; as a result, the randomness in injury occurrence does not change with time in each loading; in other words, if a dose of level x causes no injury at time t₁, a dose of the same level at a later time won’t result in injury; a static load (dose) of level x over a long time will have exactly the same effect as the same load over a short time.

A prominent property of memoryless loading process is that the injury probability is determined solely by the peak value of the loading profile: increasing the load (the dose) to level x and then releasing the load has the same effect as a static load of level x in causing injury.

A quasi-static loading test (i.e. dynamic loading with small loading rate) can be viewed approximately as a memoryless process. We compare two inference formulations for constructing the dose injury relation from the observed samples of fracture tolerance in a quasi-static loading test.

4.2. Two Inference Formulations

Method 1: Estimating $p\left(x=Y_j\right)$ and fitting logistic function to (Y_j, p(Y_j)).

This is the methodology followed by [1]. First, the measured samples of fracture tolerance are sorted into an increasing sequence {Y_j}. Under the assumption of a memoryless process, the sorted data sequence indicates that if the dose is maintained at x = Y_j, then j samples out of the m samples tested are injured (i.e., fractured). These are samples 1 through j in the sorted sequence, with observed fracture tolerance less than or equal to Y_j. Thus, at input dose x = Y_j, it is reasonable to estimate the corresponding injury probability as:

$p\left(x=Y_j\right)\approx \frac{j}{m}$

We use the logistic function form $f_L\left(x;D_{50},W\right)$ defined in (1). The dose-injury relation is described by shape parameters (D₅₀, W), which are determined by fitting function $f_L\left(x;D_{50},W\right)$ to data $\left\{\left(Y_j,j/m\right),j=1,2,\cdots ,m\right\}$.

$\text{Method1:}\left(D_{50}^{\left(est\right)},W^{\left(est\right)}\right)=\underset{\left(D_{50},W\right)}{\arg \min }{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\left(f_L\left(Y_j;D_{50},W\right)-\frac{j}{m}\right)}^2}$

This is a non-linear least squares fitting problem.

Method 2: Using the framework of dose propagation uncertainty.

With the injury function given in (2), the task of inferring the dose injury relation is reduced to estimating parameters (D₅₀, W) from measured samples of fracture tolerance: {Y_j}. In a loading process, the input dose x(t) and the target dose Z(t) are functions of time. As discussed above in deriving the dose propagation uncertainty formulation, the difference (Z(t) – Z^(c)) has the expression:

$Z\left(t\right)-Z^{\left(c\right)}=x\left(t\right)-\mu -s\xi ,\xi ~N\left(0,1\right)$

The uncertainty in dose propagation is attributed to biovariability of the test sample (i.e., bone density, bone geometry) and attributed to realization of the experimental setup (i.e., where the impactor hits and where the measurement gauge is attached). These factors are invariant with respect to time during each loading. As a result, the dose propagation uncertainty, (–μ – sξ), is independent of time in each loading. Quantity (–μ – sξ), however, is still a random variable, fluctuating among individual realizations of loading and among individual test samples.

In a memoryless process, fracture occurs when the target dose Z(t) reaches the critical threshold Z^(c). In the three-point bending tests, the input dose x(t) increases with time until fracture. Let t_f be the time of fracture. At fracture, the load and the dose propagation uncertainty are related by

$0={\left(Z\left(t\right)-Z^{\left(c\right)}\right)|}_{t=t_f}={\left(x\left(t\right)-\mu -s\xi \right)|}_{t=t_f}$

By definition, fracture tolerance Y is the observed value of x at fracture: Y = x(t_f). It follows that the fracture tolerance Y has the distribution:

$Y=x\left(t_f\right)=\mu +s\xi ~N\left(\mu ,s^2\right)$

This is the distribution of Y among individual realizations of loading and among individual test samples. The distribution of Y tells us that parameters (μ, s) are the mean and standard deviation of Y. Given data set $\left\{Y_j,j=1,2,\cdots ,m\right\}$, parameters (μ, s) are estimated using the sample mean and sample standard deviation:

$\begin{array}{l}{\mu }^{\left(est\right)}=\frac{1}{m}{\displaystyle \underset{j=1}{\overset{m}{\sum }}{Y}_j}\\ s^{\left(est\right)}=\sqrt{\frac{1}{m-1}{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\left(Y_j-{\mu }^{\left(est\right)}\right)}^2}}\end{array}$

Using the relation between the uncertainty parameters (μ, s) and the shape parameters (D₅₀, W) derived in (3), we estimate (D₅₀, W) as

$\text{Method2:}\{\begin{array}{l}{D}_{50}^{\left(est\right)}=\frac{1}{m}{\displaystyle \underset{j=1}{\overset{m}{\sum }}{Y}_j}\\ W^{\left(est\right)}=2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right)\sqrt{\frac{1}{m-1}{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\left(Y_j-{\mu }^{\left(est\right)}\right)}^2}}\end{array}$

In contrast to the non-linear least squares fitting in Method 1, Method 2 requires only straightforward calculation of sample mean and sample standard deviation. In addition, Method 2 produces the maximum likelihood estimates of parameter (D₅₀, W) under the assumption that observed samples of fracture tolerance {Y_j} follow a normal distribution.

4.3. Comparison of the Two Inference Methods

Table 1 below compares the performance of the two methods in estimating the shape parameters (D₅₀, W) from data. The data sets used in comparison study are generated by simulations using parameter values D₅₀ = 16 and W = 5. Note that mathematically we can shift and scale the dose-injury function to make D₅₀ = 0 and W = 1. The error study is not affected by the particular choice of (D₅₀, W). Each data set contains m = 10 samples of fracture tolerance, reflecting the small sample size in real bending tests. We examine the errors in the estimated median injury dose and in the estimated width:

$\left(D_{50}^{\left(est\right)}-D_{50}\right)$ and $\left(W^{\left(est\right)}-W\right)$.

Table 1. Comparison of the two inference formulations.

We calculate the mean, the standard deviation (std) and the RMS (root-mean-square) of each error using Monte Carlo simulations. The values reported in Table 1 are based on N = 500,000 independent data sets, each data set yielding an independent sample of $\left(D_{50}^{\left(est\right)},W^{\left(est\right)}\right)$. Based on the results in Table 1, we conclude that Method 2 is more accurate (having smaller error) than Method 1 in every aspect. In particular, Method 2 is unbiased for estimating the median injury dose. In addition, Method 2 is much simpler to implement computationally than Method 1.

5. Scaling of Fracture Tolerance vs Loading Rate—The Case of an Injury Process with Memory

In the previous section, we reviewed the formulation of dose propagation uncertainty and extended it to the case of quasi-static dynamic loading experiments, which are approximately memoryless. In a memoryless process, the effect of input dose is instantaneous with no delay in time, and if it does not result in injury immediately there is no partial injury or left-over effect for a later time. As a result, in a memoryless injury process, the observed fracture tolerance is not affected by the loading rate.

In order to accommodate experiments with high loading rates, in particular, in order to compare test data obtained with different loading rates, we need to relax the memoryless assumption and consider the memory effect. We study the dependence of observed fracture tolerance on loading rate. For simplicity, we neglect the randomness of fracture tolerance among individual samples, and focus solely on the effect of loading rate. In general, bone is visco-elastic-plastic, which has memory. Below, in the discussion of memory effect, we use the stress σ as the input dose. The situation of using other predictor variables as the input dose may have similar behaviors.

5.1. Exponential Model of Fracture Time vs Applied Static Stress

Under a constant applied stress σ, a visco-elastic-plastic material undergoes plastic deformation at a pace that is dependent on the applied stress σ. Given sufficiently long time, it will eventually fracture. In [7], one model for fracture is that the fracture time t_f(σ) decreases exponentially with respect to the applied static stress σ.

$t_f\left(\sigma \right)=t_0\exp \left(\frac{-\sigma }{{\sigma }_0}\right)$ (4)

where

t₀: fracture time at zero applied stress (very very large);

σ₀: stress increment that reduces the fracture time by a factor of e = 2.71828 …

We view the applied stress σ as contributing uniformly over time interval [0, t_f] in driving the plastic deformation toward the eventual fracture at t_f. We measure the contribution to fracture in percentage so that fracture occurs when the total contribution reaches 1 = 100%. With the exponential model of fracture time vs applied static stress, the contribution of stress σ in time interval $\left[t,t+\Delta t\right]$ is

$\text{contribution of}\sigma \text{in}\left[t,t+\text{d}t\right]=\frac{1}{t_f\left(\sigma \right)}\text{d}t=\frac{1}{t_0}\exp \left(\frac{\sigma }{{\sigma }_0}\right)\text{d}t$

In dynamic loading experiments, the applied stress varies with time: σ(t). The cumulative contribution of applied stress profile {σ(t)} over time interval [0, T] is

$\text{contribution of}\left\{\sigma \left(t\right),0<t<T\right\}=\frac{1}{t_0}{\displaystyle {\int }_0^T\exp \left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)\text{d}t}$

In a dynamic loading experiment designed to ensure fracture occurring within a reasonably observable time, the applied stress σ(t) increases monotonically with time t until fracture. Fracture occurs when the cumulative contribution reaches 100%. To illustrate the effect of loading rate in a simple setting, we first consider the case of σ(t) increasing linearly with time t until fracture:

$\sigma \left(t\right)=\eta t$ until fracture

where η is the stress increase rate. We introduce time scale t^*.

$t^{*}={\sigma }_0/\eta$ : time scale for increasing stress by σ₀

here σ₀ is the stress increment that increases the contribution toward fracture per unit time by a factor of e = 2.71828. Accordingly, in the dynamic loading, t^* is the time increment that increases the contribution toward fracture per unit time by a factor of e = 2.71828.

We calculate the cumulative contribution of {σ(t) = ηt} over time interval [0, T].

$\frac{1}{t_0}{\displaystyle {\int }_0^T\exp \left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)\text{d}t}=\frac{1}{t_0}{\displaystyle {\int }_0^T\exp \left(\frac{t}{t^{*}}\right)\text{d}t}=\frac{t^{*}}{t_0}\left(\exp \left(\frac{T}{t^{*}}\right)-1\right)$

Let T_f(η) denote the time to fracture under dynamic loading {σ(t) = ηt}. It satisfies:

$\frac{t^{*}}{t_0}\left(\exp \left(\frac{T_f\left(\eta \right)}{t^{*}}\right)-1\right)=1$

Solving for T_f(η) we obtain

$T_f\left(\eta \right)=t^{*}\ln \left(1+\frac{t_0}{t^{*}}\right)$

In all meaningful dynamic loading experiments, we have $T_f\left(\eta \right)\gg t_0$, which implies $t_0/t^{*}\gg \text{1}$. That allows us to write the fracture time T_f(η) approximately as

$T_f\left(\eta \right)\approx t^{*}\ln \left(\frac{t_0}{t^{*}}\right)$

The fracture tolerance, σ_f(η), is defined as the stress at fracture.

${\sigma }_f\left(\eta \right)\equiv {\sigma \left(t\right)|}_{t=T_f\left(\eta \right)}={\sigma }_0\ln \left(\frac{t_0\eta }{{\sigma }_0}\right)$

Let $\stackrel{\dot{}}{\epsilon }$ denote the strain increase rate (strain rate). Suppose the stress increase rate η is proportional to the strain rate $\stackrel{\dot{}}{\epsilon }$ until fracture, $\eta =c\stackrel{\dot{}}{\epsilon }$ until fracture.

We write fracture tolerance σ_f as a function of strain rate $\stackrel{\dot{}}{\epsilon }$

$\begin{array}{c}{\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)=\ln \left(10\right){\sigma }_0{\log }_{10}\left(\frac{t_{0}c}{{\sigma }_0}\right)+\ln \left(10\right){\sigma }_0{\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)\\ =C_0+C_1{\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)\end{array}$ (5)

With the exponential model, the theoretically predicted fracture tolerance ${\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)$ is a linear function of the logarithm of strain rate. In light of the linear relation of ${\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)$ vs ${\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)$ described in (5), we propose a loading rate adjusted dose for predicting fracture risk.

$x_{\begin{array}{c}\text{Loadingrate}\\ \text{adjusteddose}\end{array}}\equiv {\sigma }_{\text{max}}\left(\stackrel{\dot{}}{\epsilon }\right)-\ln \left(10\right){\sigma }_0{\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)$ (6)

where ${\sigma }_{\text{max}}\left(\stackrel{\dot{}}{\epsilon }\right)$ is the maximum stress over a pre-selected time period [0, T] or over the entire time course of a loading. The effect of this loading rate adjusted dose in causing fracture is expected to be independent of the loading rate when the underlying exponential model is a valid description of the fracture process. In other words, the fracture probability is solely determined by the value of the loading rate adjusted dose regardless of the loading rate used in tests. For this reason, the loading rate adjusted dose serves as a formulation for unifying test results from experiments with different loading rates.

5.2. Local Linear Relation between Stress Rate and Strain Rate

In the mathematical derivation above, we used the linear relation

$\frac{\text{d}\sigma \left(t\right)}{\text{d}t}=c\stackrel{\dot{}}{\epsilon }$ (7)

The theoretically predicted linear dependence of ${\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)$ on ${\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)$, shown in Equation (5), does not require linear relation (7) be globally satisfied during the entire loading process. To derive Equation (5), we only need linear relation (7) in a time interval near fracture that contains the dominant contribution toward fracture. The detailed analysis for this general case is presented in Appendix A.

5.3. Power Law Model of Fracture Time vs Applied Static Stress

In another model for fracture [7], the time to fracture t_f(σ) under a constant applied stress σ follows a power law

$t_f\left(\sigma \right)=t_0{\left(\frac{\sigma }{{\sigma }_0}\right)}^{-b}$ (8)

With the power law model, we show that the logarithm of fracture tolerance is a linear function of the logarithm of strain rate.

$\begin{array}{c}{\log }_{10}{\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)={\log }_{10}{\sigma }_0+\frac{1}{b+1}{\log }_{10}\left(\frac{c\left(b+1\right)t_0}{{\sigma }_0}\right)+\frac{1}{b+1}{\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)\\ =C_0+C_1{\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)\end{array}$ (9)

The derivation of relation (9) is presented in Appendix B. After the completion of the manuscript, the authors became aware that relation (9) was derived in [8] for the case where the stress rate, the strain rate and their ratio are all constants, independent of time. In Appendix B we present a derivation for a slightly more general case where the stress rate and the strain rate only need to be locally proportional in a time interval leading to fracture. Both derivations were based on the view that cumulative damage increases at a rate that is solely determined by the applied stress and is independent of damage status. Fracture occurs when the cumulative damage reaches 100%.

For the power law model, the linear relation is between ${\log }_{10}\left({\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)\right)$ and ${\log }_{10}\left(\stackrel{\dot{}}{\epsilon }\right)$. Accordingly, the loading rate adjusted dose for predicting fracture probability is

$x_{\begin{array}{c}\text{Loadingrate}\\ \text{adjusteddose}\end{array}}\equiv {\log }_{10}\left({\sigma }_{\max }\left(\stackrel{\dot{}}{\epsilon }\right)\right)-\frac{1}{b+1}{\log }_{10}(\epsilon \dot{})$

This loading rate adjusted dose is based on the power law model.

5.4. Comparison of the Exponential Model and the Power Law Model

We examine the two theoretical models above for fracture tolerance vs strain rate using the data of compression tests of human femoral cortical bone samples [6]. When compressed, respectively, along the longitudinal and the transverse directions, the fracture tolerance of the bone material varies with the applied strain rate $\stackrel{\dot{}}{\epsilon }$. Table 2 summarizes the test results reported in [6] for strain rates ranging from $\stackrel{\dot{}}{\epsilon }$ = 10^–3/s to $\stackrel{\dot{}}{\epsilon }$ = 10³/s. The results include fracture stress (i.e., stress at fracture) and fracture strain (i.e., strain at fracture).

We fit the exponential model (5) and the power law model (9) to the data of compression fracture stress vs strain rate. Figure 1 plots the data and the fittings of the two models in each of the longitudinal and transverse compressions. Figure 1 demonstrates that both models fit fairly well with the data. The two models are close to each other in the range of data. Coefficient ln(10)σ₀ in exponential model (5) and coefficient 1/(1 + b) in power law model (9) are both expected to be positive based on their physical meaning. As a result, both models predict the fracture stress as an increasing function of strain rate. Indeed, experimental results from [6] confirm this theoretical prediction. Given the fracture stress measured at a particular strain rate, we like to theoretically predict the fracture stress at a different strain rate for samples of the same type. We calculate the prediction using either the exponential model (5) or the power law model (9).

Table 2. Compression test results of human femoral cortical bone in [6].

Figure 1. Measured fracture stress vs strain rate from [6] and fittings of the two models, in each of the two compression directions: longitudinal and transverse.

Exponential model:

${{\sigma }_f|}_{\stackrel{\dot{}}{\epsilon }=r_2}={{\sigma }_f|}_{\stackrel{\dot{}}{\epsilon }=r_1}+\ln \left(10\right){\sigma }_0{\log }_{10}\left(\frac{r_2}{r_1}\right)$ (10)

Power law model:

${\log }_{10}{{\sigma }_f|}_{\stackrel{\dot{}}{\epsilon }=r_2}={\log }_{10}{{\sigma }_f|}_{\stackrel{\dot{}}{\epsilon }=r_1}+\frac{1}{b+1}{\log }_{10}\left(\frac{r_2}{r_1}\right)$ (11)

For the purpose of predicting the fracture stress based on an observed fracture stress at a particular strain rate, we only need coefficient ln(10)σ₀ in (5), or coefficient 1/(1 + b) in (9). For each of the compression directions, the values of ln(10)σ₀ in (5) and 1/(1 + b) in (9) are estimated by fitting the two models to test data. The fitting results are reported in Table 3. In Table 3, coefficient ln(10)σ₀ for the longitudinal compression is significantly higher than that for the transverse compression (27.82 MPa vs 15.38 MPa). In contrast, coefficient 1/(1 + b) in the power law model is virtually the same for the two compression directions (≈0.053), which corresponds to b ≈ 18. It will be interesting to see if exponent b in the power law model remains approximately the same for bone samples of other types besides femoral cortical bones. If that is approximately true, then the power law model is universally applicable for scaling the fracture stress vs strain

Table 3. Slope coefficients of the two models, in two compression directions.

rate regardless of bone type and regardless of compression direction. In addition to possessing a universal exponent b, the predicted fracture stress (11) based on the power law model appears mathematically more reasonable than (10) based on the exponential model. In (11), the power law model always predicts a positive value for the fracture stress σ_f while in (10), the exponential model may produce a negative value for σ_f at very low strain rate. Finally, model (11) is mathematically more appealing than model (10) in that coefficient 1/(1 + b) and all terms in (11) are dimensionless.

The derivation of exponential model (5) and power law model (9) above is based on the two corresponding models for the fracture time under a constant applied stress. The key assumption is to view the applied static stress as contributing toward fracture uniformly over the time duration up to fracture. This derivation cannot be directly extended to predict the relation between the fracture strain and the strain rate. Nevertheless, we apply the function forms of models (5) and (9) phenomenologically to fit the test results of fracture strain vs strain rate from [6]. Figure 2 plots the data and the fittings of the two models in each of the longitudinal and transverse directions.

Examining the test results from [6] and the fitting results shown in Figure 1 and Figure 2, we make several observations:

· The measured fracture stress vs strain rate is well fit by both the exponential model and the power law model in the transverse compression and in the longitudinal compression. The power law model gives a slightly better fit and yields a unified exponent b ≈ 18 in both the cases of transverse compression and longitudinal compression. This makes the power law model more universally applicable.

· The variation in measured fracture strain, as represented by error bars, is significantly larger than that of corresponding measured fracture stress. Large variation in measured fracture strain can also be expected in other types of experiments, especially in three-point bending tests. There are several reasons: 1) the spot of maximum weakness due to microstructure in a bone is at an uncertain position relative to the location of strain measurement; 2) bone inhomogeneity makes local strain strongly position dependent, varying significantly from one position to another; in contrast, the stress calculated from the applied force is fairly immune from the effect of inhomogeneity; and 3) strain measurement is affected by the uncertainty in strain gauge attachment. We will discuss more on microstructure, inhomogeneity, and uncertainty in strain gauge attachment.

· The measured fracture strain vs strain rate in transverse compression does

Figure 2. Measured fracture strain vs strain rate from [6] and fittings of the two models, in each of the two compression directions: longitudinal and transverse.

not follow any simple trend. The measured fracture strain in longitudinal compression appears to be roughly invariant with respect to the strain rate. More test results at strain rates in between 10^–3, 1, and 10³ are needed to resolve this trend. Also more tests are needed to produce reliable values with reasonably small error bars.

· Tensile test results for the effect of strain rate/dynamic loading are needed, which are probably more relevant (than the compression test results) for predicting fracture caused by bending.

6. An Injury Model Based on Existing Data of Humerus Fracture

We examine the bending test results of humerus fracture and forearm fracture from [4]. We start with the humerus data in [4]. The upper arm has only one long cylinder-shaped bone: the humerus. The measured bending moment of the upper arm complex is assumed to closely approximate the bending moment of the humerus. This allowed the authors of [4] to calculate the bending stress on the humerus using the measured bending moment of the upper arm and the bone cross-sectional properties obtained in pre-test CT scans of individual samples. The values of calculated stress vs measured strain during the initial linear response regime of the loading were used to estimate the elastic modulus [4]. The calculation of stress can be carried out all the way up to fracture. The calculated stress immediately before fracture can be used as an indirect measurement of fracture stress. In [9], this method was used to calculate the fracture stress for human femurs under impact loading. The fracture stress for humeri, however, was not calculated/reported in [4].

In contrast to the upper arm, the forearm has two bones: radius and ulna. The bending moment can only be measured for the whole forearm complex containing both bones. The bending moment of each individual bone is not directly measurable in the three-point bending tests. As a result, the bending stress of each bone cannot be determined from the available measurements.

6.1. Log-Normal Distribution of Measured Fracture Tolerance

In the bending test of [4], each humerus specimen is impacted in the direction from posterior to anterior, and the time series of bending moment and time series of strains at two locations on the humerus were recorded. The two measured strains on the humerus are:

· anterior strain > 0: tensile strain on the anterior side

· posterior strain < 0: compression strain on the posterior side

For mathematical convenience, we use the absolute value when studying the posterior strain. Let Y be the random variable modeling the measured fracture tolerance, which is defined as the peak value before fracture in a recorded time series. Here Y may be any of the 3 measured fracture tolerances: anterior strain, posterior strain or bending moment. We first find whether Y is closer to a normal distribution or to a log-normal distribution. We introduce quantity Q(Y) to distinguish these two distributions.

$Q\left(Y\right)\equiv 2{\left(\frac{E\left(Y\right)}{\text{std}\left(Y\right)}\right)}^2\left[\text{med}\left(\ln \left(Y\right)\right)-E\left(\ln \left(Y\right)\right)\right]$

where E( ) denotes the mean and med( ) the median. In Appendix C, we show that Q(Y) = 1 if Y is normal; Q(Y) = 0 if Y is log-normal. Note that a log-normal distribution has a heavy tail. If the tail of Y is heavier than that of a log-normal, then we have Q(Y) < 0. Table 4 lists the Q-values for the 3 measured fracture tolerances: anterior strain, posterior strain, and bending moment. Each Q-value is based on 11 measured samples from [4].

Examining the values reported in Table 4, we make several observations:

· Measured fracture bending moment is close to a log-normal distribution.

· For Y = measured fracture anterior strain or posterior strain, between the choices of normal and log-normal, it is more appropriate to model Y as a log-normal than as a normal. The negative value of Q indicates that the tail of Y is heavier than that of a log-normal. That is, log(Y) has a heavy tail. In other words, the fluctuations of log(strain) caused by measurement uncertainty has a heavy tail, instead of being symmetric normal.

· The bottom row in Table 4 indicates higher relative variation in the measured fracture strains than in the measured fracture bending moment, especially in the tensile anterior strain. This is consistent with the observations in split-Hopkinson pressure bar compression tests of femoral cortical bone samples [6]. Some part of the variation is due to the biovariability of samples tested, which affects all 3 measured fracture tolerances. Some other part of the variation is due to the random inhomogeneity, which makes local strain strongly position dependent and thus mainly affects measured strains. In addition, measured strains are affected by the uncertainty in measurement gauge attachment (position/orientation, distance from the bone surface). Suppose we hypothetically “homogenize” a quantity by filtering out the variations caused by random inhomogeneity and measurement uncertainty. This hypothetical

Table 4. Q-values of the 3 measured fracture tolerances in humerus data.

homogenization does not filter out the effects of biovariability in macroscopic properties (for example, variations in bone density or in bone geometry). Since the measured bending moment is fairly immune from random inhomogeneity and measurement uncertainty, it will not be significantly affected by this hypothetical homogenization. In Table 4, the measured fracture bending moment has smaller variation than the 2 measured fracture strain quantities, and is close to log-normal. Thus, it is reasonable to expect that the “homogenized” fracture strain will have smaller variation and follow the log-normal distribution. This issue will be further explored in the fluctuation/residual analysis below.

Before we end this subsection, we point out that in terms of approximate CDF, log-normal and normal distributions are not mutually exclusive. The two distributions may be practically close to each other. A random variable may be both close to a log-normal and close to a normal. This occurs when log(Y) is a normal random variable with a small standard deviation. Let $\log \left(Y\right)=\mu +\sigma \xi$ where $\xi ~N\left(0,1\right)$ and $\sigma \ll \text{1}$. We write Y as

$Y=\exp \left(\mu +\sigma \xi \right)=\exp \left(\mu \right)\exp \left(\sigma \xi \right)=\exp \left(\mu \right)\left(1+\sigma \xi +O\left({\sigma }^2\right)\right)$

which is approximately a normal random variable.

Figure 3 compares the quantile-quantile plot of bending moment vs theoretical normal and that of bending moment vs theoretical log-normal. The results

Figure 3. Quantile-quantile plot of bending moment vs theoretical normal and that of bending moment vs theoretical log-normal.

indicate that it is justified to model the fracture bending moment phenomenologically either as a normal or as a log-normal. In this study, we choose to model measured fracture tolerances as log-normal.

6.2. Fluctuations in Fracture Strain and Fracture Bending Moment of Humerus Data

Let X = anterior strain or posterior strain or bending moment. The value of X at fracture is the fracture tolerance Y. Since we model ln(Y) approximately a normal random variable, it is more appropriate to use ln(X) as the input dose in the dose propagation uncertainty formulation, instead of using X as the input dose. As we derived in Section 3, the shape parameters (D₅₀, W) in the injury model for input dose ln(X) are expressed in terms of E(ln(Y)) and std(ln(Y)), the mean and standard deviation of ln(fracture tolerance). Based on observed samples of ln(Y), we approximate its mean and standard deviation as

$\begin{array}{l}\overline{\ln \left(Y\right)}=\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\ln {\left(Y\right)}_j}\\ \text{std}\left(\ln \left(Y\right)\right)\approx \sqrt{\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left(\ln {\left(Y\right)}_j-\overline{\ln \left(Y\right)}\right)}^2}}\end{array}$

In the above, the standard deviation uses N as the denominator, instead of (N – 1). This is to make it having the same form as the RMSD (root-mean-square deviation) described below so that we can directly compare these two quantities with each other.

In [4], two attributes were given for each sample tested: age (A) and body mass (M) of the cadaver donor. We explore if some of the variation in ln(Y) can be explained by variations in age (A) and/or in body mass (M). We consider the data set

$\left\{\left(A_j,M_j,Y_j\right),j=1,\cdots ,N\right\}$

We use a linear function of (A, M) as the fitting function for ln(Y):

$f\left(A,M\right)=\left(c_0+c_{A}A+c_{M}M\right)$

We fit function f(A, M) to {ln(Y)_j}. The RMSD of the fitting is

$\text{RMSD}=\sqrt{\frac{1}{N}\underset{\left(c_0,c_A,c_M\right)}{\min }{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left(c_0+c_A{A}_j+c_M{M}_j-\ln {\left(Y\right)}_j\right)}^2}}$

The approximation of std(ln(Y)) given above is simply the RMSD of fitting using constant function $f\left(A,M\right)=c_0$. Figure 4 shows data of ln(measured fracture bending moment) and the least squares fitting of $\left(c_0+c_{A}A+c_{M}M\right)$, which yields

Fitting results for bending moment:

$c_0=4.57$,

$c_A=-2.77\times {10}^{-3}$ (p-value, 0.41),

$c_M=9.78\times {10}^{-3}$ (p-value, 0.063).

Figure 4. Data of ln(measured fracture bending moment) and least squares fitting of $\left(c_0+c_{A}A+c_{M}M\right)$, which yields $c_0=4.57$, $c_A=-2.77\times {10}^{-3}$ and $c_M=9.78\times {10}^{-3}$. Left panel: data vs age (A). Right panel: data vs body mass (M).

The negative value of coefficient c_A suggests that the fracture bending moment decreases with age while the positive value of c_M indicates that the fracture bending moment increases with body mass. The data and fitting for ln(measured fracture extension strain) are displayed in Figure 5. The least squares fitting $\left(c_0+c_{A}A+c_{M}M\right)$ for ln(strain) gives.

Fitting results for extension strain:

$c_0=-0.212$,

$c_A=1.2\times {10}^{-2}$ (p-value, 0.35),

$c_M=-2.53\times {10}^{-5}$ (p-value, 0.999).

The coefficients obtained in fitting seem to indicate that the fracture strain increases with age and is independent of body mass. Given the p-value of 0.35, the predicted trend of fracture strain increasing with age is not statistically significant. It may be an artifact of small sample size and outlier fluctuations in biovariability not correlated with age. Here we treat this fitting result only as a placeholder in the framework.

Next we study the fluctuations in ln(Y) where Y is any of the 3 measured fracture quantities. We like to examine the part of variations in ln(Y) that is not correlated with (A, M). For that purpose, we first model the dependence of ln(Y) on (A, M) using a least square fitting of $\left(c_0+c_{A}A+c_{M}M\right)$. The variation in the residual $\left(\ln \left(Y\right)-\left(c_0+c_{A}A+c_{M}M\right)\right)$ serves our purpose and is represented by the RMSD (root-mean-square deviation) of the fitting. The RMSD contains the effects of 3 factors listed in Table 5.

The effects of factors ii) and iii) vary significantly in magnitude among the 3 measured fracture quantities. Random local inhomogeneity of bone significantly increases the variation in strain while having minimal effect on stress; the location and the method of attaching strain gauges indirectly on the bone introduces additional uncertainty in measuring bone strain. In contrast, the variation in macroscopic properties is expected to have approximately the same effect on all 3 measured fracture quantities.

Figure 5. Data of ln(measured fracture extension strain) and least squares fitting of $\left(c_0+c_{A}A+c_{M}M\right)$, which yields coefficients $c_0=-0.212$, $c_A=1.2\times {10}^{-2}$ and $c_M=-2.53\times {10}^{-5}$. Left panel: data vs age(A). Right panel: data vs body mass(M).

Table 5. Three factors that contribute to RMSD of fitting.

By comparing the RMSDs of fitting for the 3 fracture quantities, we identity which quantity is least affected by factors ii) and iii) above. In this way, we identify which quantity is the best candidate for isolating and quantifying the effect of factor i).

Table 6 displays the standard deviation, the RMSD of fitting, and the predicted mean value of ln(Y) at (A = 55 y, M = 55 kg) for the 3 measured fracture quantities. In Table 6, of the 3 fracture quantities, the fracture bending moment has the smallest RMSD, further supporting the assertion that bending moment is the most accurately measured and the least affected by random local inhomogeneity and measurement uncertainty. The RMSDs of fitting for the 2 strains indicate that large fluctuations remain after the dependence on age and body mass is discounted. These large fluctuations in the two measured strains likely contain significant contributions from factors ii) and iii), in addition to the contribution from factor i). Both the random local inhomogeneity and the measurement uncertainty increase the variation in measured strain. In contrast, the fracture bending moment is calculated based on: 1) the peak impact force exerted on the bone; and 2) the bone geometry. Both 1) and 2) are fairly immune from the effects of measurement uncertainty and random local inhomogeneity.

6.3. Elasticity-Homogenized Strain in Experiments and in ATBM

Any model designed to simulate a generic (as opposed to specific) individual usually does not account for random local inhomogeneity in bone microstructure. In particular, ATBM simulations treat each bone as an object of uniform

Table 6. Comparison of fluctuations in 3 measured fracture quantities in humerus data.

material properties. Therefore, the strains calculated in ATBM finite element model (FEM) simulations are incompatible with the measured peak strains from [4]. In the ATBM-FEM of the arm, the cortical bone of the humerus is modeled as a shell of finite thickness surrounding an interior of trabecular bone and being surrounded by flesh outside. Each bone/tissue type has homogeneous material properties within its region. In the FEM, the geometry of each bone type/tissue type is accurately represented in the finite element model. The material properties, however, are set as uniform and homogeneous inside each bone type region. In a simulation, the calculated strain is the strain in a 3-D region that models the cortical bone with a realistic geometry, but is homogenized to uniform material properties. As a result, the calculated strain does not contain the effect of position-dependent fluctuations caused by random inhomogeneity in bone.

In addition, three-point bending experiments have measurement uncertainty in measured strains because they are measured at the tissue surface of the upper arm complex with strain gauge location and attachment uncertainty. This measurement uncertainty is not reflected in the FEM. Finally, the FEM simulates each bone type as perfectly elastic throughout the entire loading process, disregarding the event of fracture and the material property changes leading to fracture. In real bones, plasticity will inevitably become more prominent when the strain is large especially in the phase immediately preceding fracture. The effect of plasticity is not included in ATBM-FEM. In summary, the ATBM-FEM excludes the effects of factors ii) and iii) listed in Table 5 (i.e., random local inhomogeneity and measurement uncertainty). In addition, it excludes the effect of bone plasticity, which plays a non-negligible role in the mechanism of bone fracture. In the ATBM framework for injury risk assessment, an injury model is built on experimental injury data and ATBM simulation results are compared to the injury curve to predict injury probability. We note that existing empirical injury models were based on a straightforward interpretation of experimental injury data. They contain all of the uncertainty inherent in the data set used for building the logistic regression. But they lack the ability to parse the various contributions.

We are not aiming at predicting fracture strains in future three-point bending tests of PMHS samples where strains are measured in a similar fashion as in [4]. Rather, our goal is to integrate ATBM simulations with a properly built injury model to predict bone fracture risk caused by a blunt impact. For the purpose of integrating idealized ATBM simulations and real test data, we require additional insight into distinguishing different types of variations and dissecting their effects on the measured quantities in experiments. The ultimate goal of JNLWD’s efforts is to assess the overall risk associated with a given projectile in use against the (younger, living) population of potential targets. The specific goal of our study here is to reconstruct the injury model based on existing test data under the constraint that the injury model is compatible with the strains calculated in ATBM simulations. The injury model should also possess the capability of accommodating additional uncertainty in the input dose, in the dose propagation and in the critical threshold for injury. This capability enables us to account for both natural variation in the target population and extrapolation from the population available for data collection (elderly, deceased) to the desired target population (younger, living). For a subject impacted by a non-lethal weapon projectile, the properly built injury model allows a meaningful risk prediction from the calculated strain in ATBM simulations.

To connect real test data and idealized ATBM simulations in a proper dose-injury model, we introduce the concept of “elasticity-homogenized strain”, which is defined as what the strain would be if the three assumptions below are hypothetically satisfied during the entire loading process:

· Each bone type remains perfectly elastic for the full range of loading;

· Each bone type has homogeneous material properties;

· There is no measurement uncertainty.

The strains calculated in the ATBM-FEM are elasticity-homogenized strains. Real test data reflect all effects of factors i), ii) and iii) in Table 5. Elasticity-homogenized strains contain only the effect of factor i). To extract the essential statistics of “elasticity-homogenized strain at fracture”, we need to exclude the effects of factors ii) and iii) (random local inhomogeneity and measurement uncertainty). As we analyzed above, the measured fracture bending moment is much less affected by these two factors in comparison with the measured fracture strains. The RMSD of fitting for the fracture bending moment represents the variation in measured fracture bending moments after the dependence on age and body mass (A, M) is discounted. Thus, the RMSD of fitting for the fracture bending moment is the best candidate available in existing data for isolating and quantifying the effect of factor i) in Table 5, which includes the variation in macroscopic properties and uncertainty in dose propagation due to biovariability that is not explained by age or body mass. The relative variation of bending moment is calculated as the RMSD of fitting for ln(bending moment). We use the calculated relative variation of bending moment to estimate the relative variation of the hypothetical elasticity-homogenized strain, which is not observable in experiments. Using this approach, we build an injury model that is based on existing test data and that is compatible with the results of idealized ATBM simulations.

6.4. A Dose-Injury Model for the Standard Group of Age = 55 y and Body Mass = 55 kg

The hypothetical elasticity-homogenized strain introduced above serves as the key connection between the three-point bending test data and the ATBM simulations. Based on the test data, we construct an injury model for the standard group of age 55 y and body mass = 55 kg. The dose-response model predicts the fracture probability from the hypothetical elasticity-homogenized strain. The injury curve is specified by two shape parameters: median injury dose D₅₀ and width W. When applying the injury model to a general group of subjects, we will update the shape parameters (D₅₀, W) in the framework of dose propagation uncertainty to incorporate the effects of age and body mass (A, M) [2]. We build the injury model for the standard group in steps as described below.

· Log-normal distribution of elasticity-homogenized strain

The test data indicates that the bending moment at fracture follows a log-normal distribution. Both the bending moment and the hypothetical elasticity-homogenized strain are mainly influenced by factor i) in Table 5. It is reasonable to expect that the elasticity-homogenized strain at fracture also follow a log-normal distribution.

· The dose quantity

We use ln(elasticity-homogenized strain) as the input dose in the injury model for two reasons: 1) we use ln( ) because the strain at fracture has a log-normal distribution; and 2) we use the strain instead of the bending moment because the strain is readily available from ATBM simulations in all situations. In contrast, the calculation of bending moment depends on the particular experimental setting of three-point bending tests [4], which excludes the option of using ln(bending moment) as the input dose in the general situation of impact by a projectile.

· Critical threshold for ln(elasticity-homogenized strain)

Let Z^(c) be the critical threshold on ln(elasticity-homogenized strain) for the standard group of age = 55 y and body mass = 55 kg. Fracture occurs when ln(elasticity-homogenized strain) ≥ Z^(c). Due to biovariability and dose propagation uncertainty, there is still randomness in threshold Z^(c) even within the standard group. For a given subject and a given experimental realization, we have

Z^(c) = ln(elasticity-homogenized strain at fracture)

The statistics of Z^(c) is not readily available from test data since the hypothetical elasticity-homogenized strain is not directly measurable in experiments. To estimate the mean and standard deviation of random Z^(c), we start with the distribution of ln(measured strain at fracture) and then filter out as much as possible the effects of age and body mass and filter out the effects of random local inhomogeneity and measurement uncertainty (factors ii) and iii) in Table 5).

· Mean of the critical threshold Z^(c)

In the least squares fitting, we predict the mean of ln(strain at fracture) for subjects of age 55 y and body mass 55 kg (Table 6). We use the predicted mean of ln(strain at fracture) to model the mean of the non-observable ln(elasticity-homogenized strain at fracture) for the standard group, which is E(Z^(c)).

$E\left(Z^{\left(c\right)}\right)=\ln \left(1.563\right)=0.447\equiv z_0$

· Standard deviation of the critical threshold Z^(c)

As we discussed above, the measured bending moment and the hypothetical elasticity-homogenized strain are both fairly immune from factors ii) and iii) listed in Table 5. The RMSD of fitting ln(bending moment at fracture) to f(A, M) discounts the dependence on (A, M) and captures the variation caused by factor i), which is the biovariability not explained by (A, M). The basic idea is that when a measurable quantity and a similar non-observable quantity are both mainly affected by factor i), we simply use the variation of the measurable one to model the variation of the non-observable one. Specifically, the standard deviation of the non-observable ln(elasticity-homogenized strain at fracture), std(Z^(c)), is approximated by the RMSD of fitting the measurable ln(bending moment at fracture) to f(A, M). The RMSD value is given in Table 6.

$\begin{array}{c}\text{std}\left(Z^{\left(c\right)}\right)=\text{RMSD of fitting for ln}\left(\text{bending moment at fracture}\right)\\ =0.0\text{98}\equiv \Delta z\end{array}$

· A normal distribution injury model

For the standard group (age = 55 y, body mass = 55 kg), we model the critical threshold Z^(c) as a normal random variable

$Z^{\left(c\right)}~z_0+\Delta zN\left(0,1\right)$

We select x = ln(elasticity-homogenized strain) as the dose quantity for the injury model. The corresponding injury probability is given by function f_N(x; D₅₀, W) defined in (2).

$p\left(x\right)=f_N\left(x;D_{50},W\right)$,

x ≡ ln(elasticity-homogenized strain).

where the shape parameters (D₅₀, W) are

$\begin{array}{l}{D}_{50}=z_0=0.447,\\ W=2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right)\Delta z=0.251\end{array}$

For conciseness, we shall refer to “elasticity-homogenized strain” simply as strain since this is the strain quantity calculated in ATBM simulations. The dose-injury model is shown in Figure 6. In the left panel of Figure 6, the injury probability is plotted as a function of input dose x = ln(strain). In the right panel, the injury probability is plotted as a function of strain. This is the injury model of humerus fracture for the standard group of age 55 y and body mass 55 kg. It predicts the injury probability from the input dose of ln(elasticity-homogenized strain), which is the output of ATBM simulations. For a general group with distribution of (age, body mass) deviating from (55 y, 55 kg), the additional uncertainty needs to be incorporated in the injury model of humerus fracture, and separately incorporated in the ATBM simulations when calculating the strains on humerus. The effects of (age, body mass) distribution will be discussed later.

Figure 6. Injury model of humerus fracture based on test data from [4], for the standard group of age = 55 y and body mass = 55 kg. Left panel: injury probability vs ln(strain). Right panel: injury probability vs strain (%). Here the simple notation “strain” refers to the hypothetical “elasticity-homogenized strain”. The injury model is formulated with ln(strain) as the dose. The median injury dose and the width parameter are D₅₀ = ln(1.564) = 0.447 and W = 0.251.

7. An Injury Model Based on Existing Data of Forearm Fracture

The forearm has a more complex structure than that of the upper arm. The forearm has two large bones: radius and ulna. While both the radius and the ulna articulate with the humerus at the elbow joint, the ulna is primarily for anchoring the forearm to the humerus and the radius is responsible for rotating the forearm around its axis (the ulna). At the wrist joint, only the radius directly articulates with the carpal bones, which connect the hand to the forearm. This radius-ulna structure enables the pronator and supinator muscles to drive the forearm-wrist-hand complex to rotate around its axis swiftly, smoothly, with wide range, and against significant opposing torque. In the supinated position when forearms are placed on one’s own lap with palms facing up, the radius and ulna are nearly parallel to each other with the radius on the thumb side (outside). In the pronated position when the forearms are placed on one’s own lap with palms facing down, the wrist end of radius is rotated with respect to the ulna and crossing over the ulna. The forearm has two sides: the side corresponding to the palm side (front) of hand is called the anterior side; the side corresponding to the dorsal side (back) of hand is called the posterior side.

We use the fluctuation analysis approach developed in the previous section to study the forearm fracture test results in [4].

7.1. Fluctuation Analysis on Forearm Data

In three-point bending tests [4], seven forearm samples in the pronated position were struck in the anterior-posterior direction, that is, the palm-back direction with reference to hand. The experiment was designed to emulate the situation of a person driving a car, holding the steering wheel with palms facing forward, and being struck by an object from the front. In the pronated position, the radius is behind the ulna when viewed in the anterior-posterior direction. Impact in the anterior-posterior direction will hit the ulna first. For each forearm sample, three peak quantities were measured: radius strain at fracture, ulna strain at fracture, and the bending moment of the whole forearm at fracture. For the forearm data, we carry out the fluctuation analysis as described in the previous section. The results are reported in Table 7. Similar to the humerus data, the measured bending moment of the whole forearm at fracture has the smallest relative variation (0.183) among 3 measured quantities. After the least squares fitting to discount the dependence on age and body mass (A, M), its relative variation drops further to 0.103. Table 7 shows that the measured radius strain and ulna strain have significantly higher relative variations, reflecting the effects of random local inhomogeneity and measurement uncertainty, which are factors ii) and iii) in Table 5. The measured bending moment, on the other hand, is fairly immune from these two factors.

7.2. A Dose-Injury Model for the Standard Group of Age = 60 y and Body Mass = 55 kg

As in the situation of humerus data, the hypothetical elasticity-homogenized strain serves as a bridge for connecting the real test data and idealized ATBM simulations. We use ln(elasticity-homogenized strain) as the dose in the injury model.

Based on the data set of measured ulna fracture strain or radius fracture strain, we use the least squares fitting procedure to predict the corresponding mean of ln(fracture strain) for the standard group of age 60 y and body mass 55 kg. Here we select age = 60 y as the standard group because all but one of the forearm specimens tested were from cadaver donors above age 59 [4]. Using this data set to predict quantities at age 55 will be an extrapolation, which we wish to avoid.

We use the predicted mean of ln(fracture strain) at (age = 60 y, body mass = 55 kg) to model the mean of ln(elasticity-homogenized strain at fracture) for the standard group, which gives E(Z^(c)). Recall that Z^(c) is the critical threshold in the dose propagation uncertainty formulation. Even within the standard group, Z^(c) is a random variable due to biovariability. The predicted E(Z^(c)), as reported in Table 7,

Table 7. Comparison of fluctuations in 3 measured fracture quantities in forearm data.

is

Ulna fracture: $E\left(Z^{\left(c\right)}\right)=\ln \left(0.330\right)=-1.109\equiv z_0$

Radius fracture: $E\left(Z^{\left(c\right)}\right)=\ln \left(1.334\right)=0.2882\equiv z_0$

Following the same approach as in the previous section, we approximate the standard deviation of the non-observable ln(elasticity-homogenized strain at fracture), std(Z^(c)), using the RMSD of fitting the measurable ln(bending moment at fracture) to f(A, M). The RMSD value is given in Table 7.

Ulna fracture: $\text{std}\left(Z^{\left(c\right)}\right)=0.103\equiv \Delta z$

Radius fracture: same as ulna fracture.

With this approach, the critical threshold Z^(c) for predicting ulna fracture has the same standard deviation as that for predicting radius fracture. The mean of Z^(c), however, is different for the ulna fracture and for the radius fracture. We model the critical threshold Z^(c) as a normal random variable. For the standard group of age 60 y and body mass 55 kg, Z^(c) has the distributions below, respectively for ulna fracture and radius fracture:

$Z^{\left(c\right)}~z_0+\Delta zN\left(0,1\right)$

Ulna fracture: $Z^{\left(c\right)}~-1.109+0.103N\left(0,1\right)$

Radius fracture: $Z^{\left(c\right)}~0.2882+0.103N\left(0,1\right)$

We select x = ln(elasticity-homogenized strain) as the input dose for the injury model. The corresponding injury probability is given by function f_N(x; D₅₀, W) defined in (2).

$p\left(x\right)=f_N\left(x;D_{50},W\right)$,

x ≡ ln(elasticity-homogenized strain),

where the shape parameters are related to (z₀, Δz) by

$D_{50}=z_0,W=2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right)\Delta z$

Ulna fracture: $D_{50}=-1.109,W=0.264$

Radius fracture: $D_{50}=0.2882,W=0.264$

Again, for conciseness, we shall refer to “elasticity-homogenized strain” simply as strain. Figure 7 displays plots of injury probability vs strain, respectively, for the radius fracture (left panel) and for the ulna fracture (right panel). These are the injury modes for the standard group of age 60 y and body mass 55 kg. For a general group, the effects of (age, body mass) distribution will be incorporated in the two injury models, and separately incorporated in the ATBM simulations when calculating the strains on ulna and on radius. The effects of (age, body mass) distribution will be discussed later.

8. An Injury Risk Assessment Framework Based on Idealized ATBM-FEM Simulations and Real Test Data

We consider the problem of assessing fracture risk of the humerus. The fracture risk of the radius and the ulna in forearm can be studied in a similar way.

Figure 7. Injury models for forearm fracture based on data from [4]. Left panel: probability of radius fracture vs radius strain (%). Right panel: probability of ulna fracture vs ulna strain (%). Here the simple notation “strain” refers to the hypothetical “elasticity-homogenized strain”. The injury model is formulated with ln(strain) as the dose. For radius fracture, the median injury dose and the width parameter are D₅₀ = 0.288 = ln(1.334) and W = 0.264; for ulna fracture D₅₀ = –1.108 = ln(0.33) and W = 0.264. In the plots, injury probability is shown as a function of strain, instead of ln(strain).

8.1. Separate Injury Models For Extension Fracture and Compression Fracture

In the humerus data, there are two strains measured at fracture: anterior strain and posterior strain. In three-point bending experiments, the humerus is impacted in the posterior-anterior direction, which results in an extension strain on the anterior side and a compression strain on the posterior side. In bending experiments of humerus, at fracture, the magnitude of anterior strain (extension) is observed to be larger than that of posterior strain (compression). In stand-alone extension tests and compression tests, however, bones are observed to have smaller fracture tolerance for extension than for compression [10]. These observations indicate that in humerus bending tests, the fracture failure is primarily caused by the strain on the anterior side exceeding the extension tolerance while the strain on the posterior side is still below the compression tolerance. Thus, in bending tests, the measured anterior fracture strain serves as a proper sample of the critical threshold for extension failure while the measured posterior strain at fracture gives only an underestimated (lower bound) sample of the critical threshold for compression failure.

Note that even within the standard group of (age = 55 y, body mass = 55 kg), both the extension strain threshold and compression strain threshold are still random variables due to biovariability. As discussed in the analysis of previous sections, the critical threshold for strain has a log-normal distribution. For mathematical convenience, we consider the critical threshold Z^(c) for ln(strain), which is the input dose of the injury model. Critical threshold Z^(c) has a normal distribution. In the corresponding injury model for extension fracture or for compression fracture, the median injury dose (D₅₀) is given by E(Z^(c)) and the width (W) is related to std(Z^(c)) by

$W=2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right)\text{std}\left(Z^{(c)}\right)$

In the dose propagation uncertainty formulation, the injury model has the form of function f_N(x; D₅₀, W) defined in (2).

Extension fracture:

$p^{\left(E\right)}\left(x^{\left(E\right)}\right)=f_N\left(x^{\left(E\right)};D_{50}^{\left(E\right)},W^{\left(E\right)}\right)$,

x^(E) ≡ ln(extension strain).

Compression fracture:

$p^{\left(C\right)}\left(x^{\left(C\right)}\right)=f_N\left(x^{\left(C\right)};D_{50}^{\left(C\right)},W^{\left(C\right)}\right)$,

x^(C) ≡ ln(compression strain).

For the standard group of (age = 55 y, body mass = 55 kg), the shape parameters of the injury model for extension fracture were estimated based on the bending test data.

Extension fracture:

$D_{50}^{\left(E\right)}=0.447,W^{\left(E\right)}=0.251$

For compression fracture, the lower bound estimate of $D_{50}^{\left(C\right)}$ from the bending test data is significantly smaller than $D_{50}^{\left(E\right)}$ for extension fracture. We think this underestimated $D_{50}^{\left(C\right)}$ is too low. In general, we expect the compression tolerance should be larger than the extension tolerance: $D_{50}^{\left(C\right)}>D_{50}^{\left(E\right)}$. As a reasonable placeholder until a relevant estimate becomes available, we shall use

$D_{50}^{\left(C\right)}=D_{50}^{\left(E\right)}=0.447,W^{\left(C\right)}=W^{\left(E\right)}=0.251$

The age and body mass affect both the shape parameters of injury models and the ATBM-FEM simulation setup. Their effects will be incorporated after we build the structure of injury assessment formulation by connecting the ATBM simulations and the empirical injury models.

8.2. A Formulation for Assessing the Injury Risk from a Given ATBM-FEM Simulation Setup

In this section, we consider a method for integrating the injury model into ATBM-FEM simulations and implementing the uncertainty propagation calculation.

To assess the bone fracture risk when a subject is hit by a blunt projectile, we start by assembling a random realization of the projectile-body-impact. Let Θ_ATBM denote the parameter set describing the projectile-body-impact. Basically, parameter set Θ_ATBM specifies the ATBM simulation setup.

Θ_ATBM = {parameters in ATBM setup}

Θ_ATBM consists of a large number of parameters, including but not limited to:

· mass, shape and material properties of the projectile;

· impact velocity;

· impact location and angle on the body;

· geometry of individual body components (bones, organs, tissues);

· material properties of individual body components.

Many parameters in Θ_ATBM are not deterministic. Rather, they are random variables in a random realization of ATBM setup. In ATBM simulations, these parameters could be sampled from distributions defined by characteristic parameters of the problem. For example, the impact location/angle on the body could be sampled from an empirical distribution; the elasticity modulus of humerus bone could be sampled from a distribution constructed based on the age and body mass of the subject. Once the parameters in set Θ_ATBM are sampled in a realization of ATBM setup, the strain/stress distribution can be computed over the spatial dimensions of each body component and over the entire time course of the hypothetical perfectly elastic impact. In ATBM simulations, body components are assumed to remain elastic throughout the loading; fracture failure is not directly reflected in the simulation. For a bone with homogeneous material properties, if fracture occurs, it will be at the location of maximum extension strain or at the location of maximum compression strain. Let

$\epsilon \left(\stackrel{\to }{u},t;{\Theta }_{\text{ATBM}}\right)=\text{strain distribution of humerus over}\left(\stackrel{\to }{u},t\right)$

where $\left(\stackrel{\to }{u},t\right)$ denotes the spatial and time coordinates. The strain distribution depends on the ATBM setup parameter set Θ_ATBM. For example, a larger impact velocity induces a larger strain on the bone while a larger body mass may mitigate the impact (due to the amount of muscle and tissue surrounding the bone). For simplicity, we denote the strain distribution concisely as $\epsilon \left(\stackrel{\to }{u},t\right)$. In the strain distribution, positive values are extension strain while negative values are compression strain. For compression strain, we take the maximum of its absolute value. Over the humerus bone and over time, the maximum extension strain ${\epsilon }_{\max }^{\left(E\right)}$ and the maximum compression strain ${\epsilon }_{\max }^{\left(C\right)}$ are, respectively,

$\begin{array}{l}{\epsilon }_{\max }^{\left(E\right)}=\underset{\left(\stackrel{\to }{u},t\right)}{\max }\left(\epsilon \left(\stackrel{\to }{u},t\right),0\right)\\ {\epsilon }_{\max }^{\left(C\right)}=\underset{\left(\stackrel{\to }{u},t\right)}{\max }\left(-\epsilon \left(\stackrel{\to }{u},t\right),0\right)\end{array}$

Note that although it is not explicitly shown in the concise notation, both ${\epsilon }_{\max }^{\left(E\right)}$ and ${\epsilon }_{\max }^{\left(C\right)}$ are functions of the ATBM setup Θ_ATBM.

We calculate the fracture risk corresponding to the strain distribution $\epsilon \left(\stackrel{\to }{u},t\right)$ as

$\begin{array}{l}p\left(\left\{\epsilon \right\}\right)=\max \left(p^{\left(E\right)},p^{\left(C\right)}\right)\\ p^{\left(E\right)}=f_N\left(x^{\left(E\right)};D_{50}^{\left(E\right)},W^{\left(E\right)}\right)\\ p^{\left(C\right)}=f_N\left(x^{\left(C\right)};D_{50}^{\left(C\right)},W^{\left(C\right)}\right)\end{array}$ (12)

where

$\begin{array}{l}{x}^{\left(E\right)}=\ln \left({\epsilon }_{\max }^{\left(E\right)}\right)=\ln \left(\underset{\left(\stackrel{\to }{u},t\right)}{\max }\left(\epsilon \left(\stackrel{\to }{u},t\right),0\right)\right),\\ x^{\left(C\right)}=\ln \left({\epsilon }_{\max }^{\left(C\right)}\right)=\ln \left(\underset{\left(\stackrel{\to }{u},t\right)}{\max }\left(-\epsilon \left(\stackrel{\to }{u},t\right),0\right)\right)\end{array}$

The mathematical approach of taking the larger of the extension fracture risk and the compression fracture risk is based on the view that the random biovariability in bone properties tends to influence the extension fracture threshold and the compression fracture threshold in a similar way. For example, a lower bone density is likely to decrease both the extension and compression fracture thresholds.

In the unlikely situation where the extension fracture threshold and the compression fracture threshold are independent of each other, the fracture risk is governed by

$1-p\left(\left\{\epsilon \right\}\right)=\left(1-p^{\left(E\right)}\right)\left(1-p^{(C)}\right)$

which leads to

$p\left(\left\{\epsilon \right\}\right)=p^{\left(E\right)}+p^{\left(C\right)}-p^{\left(E\right)}\cdot p^{(C)}$

We think it is unreasonable to assume the independence of extension and compression fracture thresholds. So, we shall use (12), which is based on the more reasonable assumption that the extension and compression fracture thresholds are highly correlated in the randomness of biovariability.

Injury model (12) calculates the fracture probability corresponding to a high-dimensional input in the form of strain distribution $\epsilon \left(\stackrel{\to }{u},t\right)$ over space and time. For each sampled realization of ATBM setup Θ_ATBM, the output of ATBM simulations gives the corresponding strain distribution $\epsilon \left(\stackrel{\to }{u},t\right)$. In this way, we have constructed a computational framework that maps a sampled realization of ATBM setup Θ_ATBM to the corresponding probability of bone fracture p.

$\begin{array}{l}\underset{\begin{array}{c}\text{ATBM}\\ \text{setup}\end{array}}{\underset{︸}{{\Theta }_{\text{ATBM}}}}\stackrel{\text{ATBMsimulation}}{\to }\underset{\begin{array}{c}\text{Strain}\\ \text{distribution}\end{array}}{\underset{︸}{\left\{\epsilon \left(\stackrel{\to }{u},t\right)\right\}}}\stackrel{\text{Calculateinputdose}}{\to }\underset{\begin{array}{c}\text{Doseforpredictingextension}\\ \text{fracture,compressionfracture}\end{array}}{\underset{︸}{\left(\ln \left({\epsilon }_{\max }^{\left(E\right)}\right),\ln \left({\epsilon }_{\max }^{\left(C\right)}\right)\right)}}\\ \stackrel{\text{Injurymodel}}{\to }\underset{\begin{array}{c}\text{Probabilityofextension}\\ \text{fracture,compressionfracture}\end{array}}{\underset{︸}{\left(p^{\left(E\right)},p^{\left(C\right)}\right)}}\stackrel{}{\to }\underset{\begin{array}{c}\text{Probability}\\ \text{offracture}\end{array}}{\underset{︸}{p}}\end{array}$

The computational framework summarized above predicts the fracture injury risk for the standard group of (age = 55 y, body mass = 55 kg) and for a given ATBM setup Θ_ATBM.

8.3. Effects of Age and Body Mass

Let A = age, M = body mass. We first look at how to incorporate the effect of (A, M) into the dose-injury relation. In the injury model, ln(elasticity homogenized strain) is selected as the input dose. In the dose propagation uncertainty formulation, the injury model is specified by the two shape parameters: the median injury dose D₅₀ and the width W. The median injury dose D₅₀ is the mean of the non-observable ln(elasticity homogenized strain) at fracture. D₅₀ is estimated by fitting measured ln(Y) to a linear function of (A, M) where Y is the anterior (extension) strain at fracture. The least squares fitting reveals the dependence of median injury dose D₅₀ on age (year) and body mass (kg).

Median injury dose vs (A, M) for humerus fracture:

$D_{50}^{\left(E\right)}\left(A,M\right)=0.447+1.2\times {10}^{–2}\left(A-55\right)-2.53\times {10}^{–5}\left(M-55\right)$

Here $D_{50}^{\left(E\right)}$ is for ln(elasticity homogenized extension strain). It needs to be pointed out that the dependence of $D_{50}^{\left(E\right)}$ on (A, M) given above is based on fitting to a very small set of test data. It should be regarded only as a placeholder until a better candidate is available. A more accurate dependence relation needs to be constructed based on more experimental measurements if available, or based on theoretical/empirical knowledge, or based on a combination of the two.

The fitting above based on limited data does not describe how the width of injury model varies with age and body mass. We shall assume that the width for the standard group (estimated in the previous sections) is valid for all (A, M).

Width parameter vs (A, M) for humerus fracture:

$W^{\left(E\right)}\left(A,M\right)=0.251$ for all (A, M)

Again, this serves only as a placeholder until a better candidate is available.

Next we incorporate the effect of age and body mass into the ATBM simulations. We write the ATBM setup parameter set Θ_ATBM explicitly as

${\Theta }_{\text{ATBM}}\left(A,M,\omega \right)$

In the function form above, (A, M) describes the distribution from which parameter set Θ_ATBM is sampled. But (A, M) does not completely specify Θ_ATBM. Even when (A, M) is fixed, Θ_ATBM is still a random variable. There is biovariability not explained by (A, M). In addition, the impact velocity, impact location and angle of the projectile have uncertainty. The overall randomness in Θ_ATBM is symbolically represented by ω.

8.4. A Framework for Assessing the Injury Risk of a Subject of Given Age and Body Mass

Consider a random subject of given age and body mass (A, M). We study the humerus fracture risk of the subject caused by the blunt impact of a projectile. Even though (A, M) is given, there are randomness in the subject’s biovariability and randomness in the projectile-body-impact. As a result, the injury probability varies among individual subjects and among individual impact realizations of a fixed subject. The average injury probability for a random subject of given (A, M) may be calculated following the steps below:

1) Sample the parameter set Θ_ATBM(ω) to set up ATBM simulations, based on the given value of (A, M) and other characteristic parameters of the problem;

2) Calculate the input dose for the injury model via ATBM simulations:

· run ATBM simulations with the sampled realization of Θ_ATBM(ω) to calculate the strain distribution over space and time $\epsilon \left(\stackrel{\to }{u},t\right)$ ;

· find the maximum extension strain and the maximum compression strain from the ATBM simulation results, and calculate the input dose for the injury model;

dose for predicting extension fracture:

$x^{\left(E\right)}=\ln \left({\epsilon }_{\max }^{(E)}\right)$

dose for predicting compression fracture:

$x^{\left(C\right)}=\ln \left({\epsilon }_{\max }^{(C)}\right)$

3) set up the injury model by calculating shape parameters (D₅₀, W) corresponding to the given (A, M);

extension fracture:

$D_{50}^{\left(E\right)}\left(A,M\right),W^{\left(E\right)}\left(A,M\right)$

compression fracture:

$D_{50}^{\left(C\right)}\left(A,M\right),W^{\left(C\right)}\left(A,M\right)$

4) apply each injury model on each input dose to calculate the extension fracture probability, the compression fracture probability, and the injury risk corresponding to the sampled realization of Θ_ATBM(ω);

$\begin{array}{l}{p}^{\left(E\right)}=f_N\left(x^{\left(E\right)};D_{50}^{\left(E\right)},W^{\left(E\right)}\right)\\ p^{\left(C\right)}=f_N\left(x^{\left(C\right)};D_{50}^{\left(C\right)},W^{\left(C\right)}\right)\\ p\left(\omega \right)=\max \left(p^{\left(E\right)},p^{(C)}\right)\end{array}$

5) repeat steps 1-4 above, and average the injury risk over independent realizations of ATBM setup Θ_ATBM(ω).

${p|}_{\left(A,M\right)}=E\left(p(\omega )\right)$

Mathematically, the risk assessment framework is concisely written as

1) $\left(A,M\right)\stackrel{\begin{array}{l}\text{Sampleparameters}\\ \text{tosetupATBM}\end{array}}{\to }{\Theta }_{\text{ATBM}}(\omega )$

2) $\begin{array}{l}{\Theta }_{\text{ATBM}}\left(\omega \right)\stackrel{\begin{array}{c}\text{ATBMsimulationstocalculate}\\ \text{straindistribution}\end{array}}{\to }\left\{\epsilon \left(\stackrel{\to }{u},t\right)\right\}\left(\omega \right)\\ \stackrel{\begin{array}{c}\text{Findmaximumstrain}\\ \text{tocalculatethedose}\end{array}}{\to }{x}^{\left(E\right)}\left(\omega \right)=\ln \left({\epsilon }_{\max }^{\left(E\right)}\right),x^{\left(C\right)}\left(\omega \right)=\ln \left({\epsilon }_{\max }^{(C)}\right)\end{array}$

3) $\left(A,M\right)\stackrel{\begin{array}{c}\text{Calculateshapeparameters}\\ \text{ofinjurymodel}\end{array}}{\to }\left(D_{50}^{\left(E\right)},W^{\left(E\right)}\right),\left(D_{50}^{\left(C\right)},W^{(C)}\right)$

4) $\left[\begin{array}{l}{p}^{\left(E\right)}=f_N\left(x^{\left(E\right)}\left(\omega \right),D_{50}^{\left(E\right)},W^{\left(E\right)}\right)\\ p^{\left(C\right)}=f_N\left(x^{\left(C\right)}\left(\omega \right),D_{50}^{\left(C\right)},W^{\left(C\right)}\right)\end{array}\right]\stackrel{\begin{array}{c}\text{Calculateinjuryrisk}\\ \text{forrealization}\omega \end{array}}{\to }p(\omega )$

5) $p\left(\omega \right)\stackrel{\begin{array}{c}\text{Averageinjuryriskover}\\ \text{independentrealizations}\end{array}}{\to }{p|}_{\left(A,M\right)}$ (13)

The computational framework summarized above predicts the fracture injury risk for a random subject of given age and body mass and for a random ATBM setup Θ_ATBM.

8.5. Options for Reducing the Number of ATBM Runs When Assessing the Injury Risk of a Population

Risk assessment calculation based on framework (13) requires a significant number of ATBM-FEM runs for a subject of given age and body mass (A, M). Each ATBM-FEM run is a time evolution simulation of the three-dimensional projectile-body complex using a finite element discretization. Computationally a single ATBM-FEM run is already substantially expensive. The calculation of average injury risk for a subject of given age and body mass is based on many runs with sampled realizations of Θ_ATBM(ω). If, in addition, we need to average the injury risk over many subjects with diverse values of (A, M), then the overall calculation of the average injury risk for a population will entail a large number of runs. Computationally this may be prohibitively expensive. There are two possible avenues for reducing the number of ATBM runs needed.

Option 1: Uncertainty propagation from ATBM setup to ATBM output

The ATBM-FEM simulation setup is specified by parameter set Θ_ATBM(ω), which is a random variable even at a fixed (A, M). For each realization of Θ_ATBM(ω), the corresponding input dose for the injury model is calculated via ATBM simulations. Let Dose(Θ_ATBM(ω)) denote the dose corresponding to ATBM-FEM setup Θ_ATBM(ω). Mathematically, this is represented by

$\begin{array}{l}{\Theta }_{\text{ATBM}}\left(\omega \right)\stackrel{\begin{array}{c}\text{ATBM}\\ \text{simulation}\end{array}}{\to }\text{Dose}\left({\Theta }_{\text{ATBM}}\left(\omega \right)\right)\\ \left|\right|\\ \left(\ln \left({\epsilon }_{\max }^{\left(E\right)}\right),\ln \left({\epsilon }_{\max }^{\left(C\right)}\right)\right)(\omega )\end{array}$

At a fixed (A, M), we use the median of Θ_ATBM(ω) as a representative ATBM-FEM setup.

${\Theta }_{\text{ATBM}}^{\left(0\right)}=\text{median}\left({\Theta }_{\text{ATBM}}(\omega )\right)$

We view the corresponding dose, $\text{Dose}\left({\Theta }_{\text{ATBM}}^{\left(0\right)}\right)$, as a representative dose of the impact situation. As the random setup Θ_ATBM(ω) fluctuates around the representative setup ${\Theta }_{\text{ATBM}}^{\left(0\right)}$, it is reasonable to expect that the corresponding fluctuations of Dose(Θ_ATBM(ω)) around $\text{Dose}\left({\Theta }_{\text{ATBM}}^{\left(0\right)}\right)$ are normally distributed.

$\text{Dose}\left({\Theta }_{\text{ATBM}}\left(\omega \right)\right)=\text{Dose}\left({\Theta }_{\text{ATBM}}^{\left(0\right)}\right)+N\left({\mu }_{\Theta },{\sigma }_{\Theta }^2\right)$

Furthermore, it is likely that the parameters $\left({\mu }_{\Theta },{\sigma }_{\Theta }^2\right)$ are approximately invariant with respect to (A, M), or they follow a simple function of (A, M).

$\left({\mu }_{\Theta },{\sigma }_{\Theta }^2\right)$ = invariant with respect to (A, M)

Suppose this conjecture is valid and we pre-calculate the value of $\left({\mu }_{\Theta },{\sigma }_{\Theta }^2\right)$ by carrying out repeated runs with randomly sampled ATBM-FEM setups. Once the value of $\left({\mu }_{\Theta },{\sigma }_{\Theta }^2\right)$ is known, at each (A, M) we only need a single ATBM run with the representative setup, ${\Theta }_{\text{ATBM}}^{\left(0\right)}|\left(A,M\right)$, to capture the distribution of Dose(Θ_ATBM(ω))|(A, M). There is no need to carry out multiple ATBM runs to generate samples of Dose(Θ_ATBM(ω))|(A, M). For the purpose of calculating the injury risk, we don’t even need to apply the dose-injury relation to individual samples of Dose(Θ_ATBM(ω))|(A, M) and then average over them. In the dose propagation uncertain formulation, uncertainty in the input dose Dose(Θ_ATBM(ω))|(A, M) is conveniently incorporated into the injury model by updating the shape parameters (D₅₀, W).

$\begin{array}{l}{D}_{50}^{\left(\text{updated}\right)}=D_{50}-{\mu }_{\Theta }\\ W^{\left(\text{updated}\right)}=\sqrt{W^2+{\left(2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right){\sigma }_{\Theta }\right)}^2}\end{array}$

With the updated shape parameters, the average injury probability over Dose (Θ_ATBM(ω)) is expressed in terms of the representative dose $\text{Dose}\left({\Theta }_{\text{ATBM}}^{\left(0\right)}\right)$ as

$E\left(p\left(\omega \right)\right)=f_N\left(\text{Dose}\left({\Theta }_{\text{ATBM}}^{\left(0\right)}\right);D_{50}^{\left(\text{updated}\right)},W^{\left(\text{updated}\right)}\right)$

After the injury probability is calculated separately for each of the extension fracture and compression fracture, the larger one of the two probabilities gives the injury probability of bone fracture.

In summary, by integrating ATBM-FEM simulations with the versatile dose propagation uncertainty formulation, the risk assessment of a subject requires only a single ATBM-FEM run to calculate the representative dose based on the representative setup ${\Theta }_{\text{ATBM}}^{\left(0\right)}$. Uncertainty in the ATBM setup Θ_ATBM(ω) leads to uncertainty in the input dose Dose(Θ_ATBM(ω)). This uncertainty is readily combined with other uncertainties and is reflected in the updated shape parameters of the injury model, using the dose propagation uncertainty formulation. There is no need to run Monte Carlo style simulations to model this uncertainty. The injury risk is calculated by applying the updated injury model on the representative dose.

Option 2: Sample Θ_ATBM(ω) and (A, M) jointly

Option 1 above reduces the computational cost to a single ATBM run for assessing the injury risk of each subject. It is based on the assumption that the dose corresponding to random ATBM setup Θ_ATBM(ω) is normally distributed around the representative dose.

$\text{Dose}\left({\Theta }_{\text{ATBM}}\left(\omega \right)\right)=\text{Dose}\left({\Theta }_{\text{ATBM}}^{\left(0\right)}\right)+\text{normaldistribution}$.

If this assumption is invalid, the risk assessment for a given subject following framework (13) will require multiple ATBM-FEM runs. For a population, however, we can still reduce the total number of ATBM runs needed in risk assessment computation.

Suppose we need to average the injury risk over many subjects with diverse values of (A, M) and over random realizations of ATBM setup Θ_ATBM(ω). A straightforward method completes the job in two steps:

· For each subject of a given (A, M), draw samples of ATBM setup Θ_ATBM(ω). For each sampled realization of Θ_ATBM(ω), calculate the corresponding injury risk p(ω)|(A, M). Then average the injury risk over independent realizations of Θ_ATBM(ω) to calculate p|(A, M), the injury risk for a subject of the given (A, M).

· Repeat the computation for each individual subject sampled from the (A, M) distribution of the population. Average the injury risk p|(A, M) over sampled subjects.

Overall, we need to take average of the injury risk, p(ω)|(A, M) with respect to ω and then with respect to individual values of (A, M) in the population. The space of (ω, A, M) is high dimensional. Nevertheless, p(ω)|(A, M) is a scalar quantity. The power of Monte Carlo integration over a high dimensional space is to sample it jointly instead of sampling each dimension or each subspace separately. Following this general rule of Monte Carlo sampling, we view individual values of (A, M) in the population as samples of a random variable. We sample (A, M, ω) jointly

$\left\{\left(A_j,M_j,{\omega }_j\right),j=1,2,\cdots ,N\right\}$

The average injury risk of the population is approximated as

$\underset{\begin{array}{c}\text{True}\\ \text{average}\end{array}}{\underset{︸}{\text{Average}p}}=\underset{\begin{array}{c}\text{MonteCarlo}\\ \text{approximation}\end{array}}{\underset{︸}{\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}p\left(A_j,M_j,{\omega }_j\right)}}}+\underset{\begin{array}{c}\text{Approximation}\\ \text{error}\end{array}}{\underset{︸}{O\left(1/\sqrt{N}\right)}}$ (14)

Note that the approximation error in Monte Carlo integration is proportional to $1/\sqrt{N}$, regardless of the number of dimensions in the space of (A, M, ω).

In summary, when we do need to use Monte Carlo simulations to average over many random factors, we do so by sampling the high-dimensional space jointly, instead of sampling each subspace separately.

8.6. A Computational Framework for Assessing the Injury Risk and Its Uncertainty for a Population

We consider a population with a given distribution of (A, M). We study the average injury risk of bone fracture per subject in the population. The average injury risk per subject, denoted by p_avg, is interpreted as follows. For a crowd, (n × p_avg) gives the expected number of injured in the crowd where n is the number of subjects hit by projectiles.

The mathematical meaning of p_avg is clear and is expressed in (14): p_avg = injury risk averaged over all subjects and over all random realizations of ATBM-FEM setup.

p_avg is calculated following the steps below.

Step 1. Sample (A, M)

Sample age and body mass (A, M) from the given distribution of (A, M).

Step 2. Sample Θ_ATBM(ω)|(A, M)

Sample the ATBM setup Θ_ATBM(ω)|(A, M) as follows:

a) Determine or sample the geometry and properties of the projectile based on the given projectile type description;

b) Determine or sample the geometry and properties of bones and tissues based on the sampled value of (A, M) from step 1;

c) Determine or sample the impact velocity based on the given projectile type and distance description (for example, a subject is hit by a projectile, of a certain type, launched from 100 - 200 meters away);

d) Determine or sample the impact location and impact angle on the subject body based on the given impact description (for example, a subject is hit on the left upper torso by a projectile of …).

Step 3. Calculate the input dose for predicting the fracture risk

Calculate the input dose for the injury model via ATBM simulations as follows:

a) Run ATBM simulations with setup ΘATBM(ω)|(A, M) to calculate the strain distribution over space and time;

b) Find the maximum extension strain and the maximum compression strain from the ATBM simulation results.

Step 4. Set up the dose-injury relation for predicting the fracture risk

Set up the injury model by calculating shape parameters D₅₀(A, M) and W(A, M) according to the sampled value of (A, M) from step 1.

Step 5. Calculate the injury risk corresponding to the sampled (A, M, ω)

Calculate the injury risk corresponding to the sampled (A, M, ω) by applying the dose-injury relation on the input dose generated by the ATBM simulations based on the sampled (A, M, ω).

Step 6. Repeat sampling (A, M, ω), and average the injury risk

Repeat steps 1-5 above, average the injury risk over independently sampled realizations of (A, M, ω), and calculate the uncertainty in the average injury risk. For example, we may report the average injury risk in the form of 6% ± 2%.

Note that the framework described above can readily accommodate the effect of gender by sampling (A, M, G) where G denotes the gender of a subject. Accordingly, in Step 2 above, we need a component model that describes the dependence of bone geometry and properties on (A, M, G); and in Step 4, we need a component model mapping (A, M, G) to shape parameters (D₅₀, W) in the dose-injury function. These component models are to be developed.

In summary, the computational framework described above assesses the injury risk of a population by 1) sampling ATBM setup; 2) running ATBM simulations to calculate elasticity-homogenized strain; 3) properly interpreting test data to formulate an injury model that maps the elasticity-homogenized strain to the bone fracture probability; and 4) applying the injury model to calculate the fracture probability and then averaging over independent samples of ATBM setup. The key for properly connecting real test data and ATBM simulations is the concept of elasticity-homogenized strain.

9. Concluding Remarks and Future Work

We have constructed a computational framework for assessing the bone fracture risk of a subject hit by a blunt impact projectile. The framework unifies real test data and ATBM simulations by introducing the hypothetical elasticity-homogenized strain as the input dose for predicting injury risk, and by formulating the dose-injury relation based on dose propagation uncertainty. The elasticity-homogenized strain corresponds to the strain quantity calculated in ATBM simulations where all bones and tissues are assumed to be perfectly elastic and have homogeneous material properties in each part. The very first step in integrating ATBM simulations and real test data is to recognize that the strain calculated in ATBM simulations is different from the strain measured in experiments of real bone samples. This is especially true when we look at the uncertainty and fluctuations in measured strains. As a result, mathematically it is improper to use the calculated strain directly in the place of the measured strain in a straightforward empirical injury model. In three-point bending tests, the fracture strain is measured along with the fracture bending moment. The measured fracture strain is highly affected by the random local inhomogeneity of bone and affected by the uncertainty in measurement gauge attachment location/orientation. On the other hand, the measured fracture bending moment is fairly immune from the effects of these two factors and consequently it has much smaller relative fluctuations in the test data. The fracture bending moment is unique to the situation of three-point bending tests, and thus, is unsuitable for serving as the input dose for predicting injury risk in a general blunt impact situation. However, the measured fracture bending moment provides valuable information about the effect of biovariability in the hypothetical absence of random local inhomogeneity and measurement uncertainty. This information enables us to filter out the effects of random local inhomogeneity and measurement uncertainty in the measured fracture strain. The filtering yields statistics of the non-observable elasticity-homogenized strain at fracture, which is the hypothetical strain at fracture if bones and tissues would remain perfectly elastic with homogeneous material properties all the way up to fracture. That is why it is called a hypothetical quantity. The dose-injury relation based on the hypothetical elasticity-homogenized strain serves as a bridge for connecting the real test data and the idealized ATBM simulations in a proper mathematical setting.

The measured strain quantities are naturally affected by random local inhomogeneity and measurement uncertainty. The key step in integrating test data with ATBM simulations is to interpret test data in terms of the hypothetical elasticity-homogenized strain and subsequently construct an empirical injury model with elasticity-homogenized strain as the input dose. The resulting computational framework performs risk assessment in a sequence of steps: 1) sample random realizations of ATBM setup; 2) carry out an ATBM run for each setup realization; 3) use the strain distribution from ATBM simulations to calculate the input dose for predicting the fracture risk; 4) apply the empirical injury model on the input dose to calculate the injury risk corresponding to each realization of ATBM setup, and finally; 5) examine the distribution of injury risk over individual realizations, calculate the average injury risk and calculate the uncertainty in the injury risk.

If we only need to estimate the average injury risk, we may reduce the computational cost dramatically by carrying out only a single ATBM run with the median representative setup. The uncertainty in the calculated dose corresponding to uncertainty in ATBM setup is included in the risk assessment formulation by updating the shape parameters (median injury dose and width) of the dose-injury relation. This computationally efficient approach is especially appropriate and practical in situations where the distribution of ATBM setup is only vaguely characterized with significant randomness, rather than precisely described.

The computational framework developed in this study has a robust theoretical structure for integrating ATBM simulations and real test data in risk assessment calculation. There are, however, many component models in the framework that need to be developed and refined. For example, in the dose-injury relation, the dependence of shape parameters on age and body mass is tentatively determined by fitting to bending test results of 10 humerus samples from cadaver donors of average age 55 (out of the 12 samples tested, only 10 had valid measurements in all aspects) [3] [4]. The general trend determined in the fitting appears intuitively reasonable: lower age and/or higher body mass lead to higher injury threshold and thus, lower injury risk at a given input dose. Nevertheless, the small sample size and the narrow and old age range make the obtained dependence inadequate if we use it to predict the median injury dose for a subject of age 20. This empirical dependence on age and body mass is only meant to be a placeholder until a more sophisticated component model is developed. Similarly, in the ATBM simulation setup, the dependence of bone geometry and material properties on age and body mass needs to be carefully investigated in a component model.

The computational framework described above integrates real test data and idealized ATBM simulation via the hypothetical elasticity-homogenized strain, which is a direct output from ATBM simulations and which serves as the input dose in the injury model for predicting bone fracture risk. The most challenging part is to filter out the effects of random local inhomogeneity and measurement uncertainty, and extract the statistics of the hypothetical elasticity-homogenized strain from measured strains in test data. In this study, we accomplished this filtering by comparing the distribution of measured strain with that of measured bending moment, which is fairly immune from random local inhomogeneity and measurement uncertainty. As a future project, we may consider adopting the stress as the input dose for predicting bone fracture risk. When a load is applied, the strain distribution is highly affected by the local inhomogeneity and the measured strain is further affected by the strain gauge attachment uncertainty. In contrast, in standard compression and tensile tests, the stress calculated from the applied force already reflects the homogenized stress, and thus, is fairly immune from the local inhomogeneity. In addition, the stress calculated from the applied force has minimal measurement uncertainty. We propose to adopt the stress as the input dose in a future version of the risk assessment framework. The key component of the framework is an injury model that maps the stress to the fracture probability. We propose to construct this injury model based on samples of fracture stress measured in standard compression and tensile tests of relevant bones. The compression tests on human femoral cortical bone conducted in [6] provide part of the data. More standard compression tests on other bone types and standard tensile tests are needed for building an adequate injury model. The collection of these existing and future standard compression/tensile tests on various bone types will greatly advance our understanding of bones under mechanical loads.

Disclaimer

The authors thank the Joint Non-Lethal Weapons Directorate of U.S. Department of Defense and the Naval Postgraduate School for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the U.S. Government.

Appendix A: Fracture Tolerance vs Strain Rate in the Exponential Model

We consider the model in which fracture time t_f(σ) decreases exponentially with the applied static stress σ:

$t_f\left(\sigma \right)=t_0\exp \left(\frac{-\sigma }{{\sigma }_0}\right)$

We study the case of dynamic loading where the applied stress increases with the time until fracture. In the main text, we derived the linear relation of fracture tolerance (stress) vs logarithm of applied strain rate when the stress rate is proportional to the strain rate. Specifically, we derived (5) based on (7). Here we show that (5) actually does not require (7) be globally satisfied during the entire loading process. Linear relation (5) is valid when the following conditions are satisfied:

· The stress vs time profile is linear in a time interval near fracture;

· This time interval contains dominant contribution toward fracture;

· Contribution toward fracture per unit time increases significantly over this time interval; and

· In this time interval, the stress increase rate is proportional to the strain rate.

Mathematically, these conditions are described as

1) σ(t) is linear for t > t₁ until fracture.

$\sigma \left(t\right)=\{\begin{array}{ll}w\left(t\right),\hfill & t<t_1\hfill \\ w\left(t_1\right)+\eta \left(t-t_1\right),\hfill & t>t_1\text{untilfracture}\hfill \end{array}$

2) The contribution toward fracture from $\left\{\sigma \left(t\right),0<t<t_{\text{1}}\right\}$ is negligible.

$\frac{1}{t_0}{\displaystyle {\int }_0^{t_1}\exp \left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)\text{d}t}\ll 1$

3) Let T_f = the fracture time. The contribution per unit time at T_f is significantly larger than that at t₁.

For t ≥ t₁, the contribution toward fracture per unit time is

$\frac{1}{t_0}\exp \left(\frac{w\left(t_1\right)+\eta \left(t-t_1\right)}{{\sigma }_0}\right)$

The ratio of the contribution at T_f and that at t₁ gives us

$\exp \left(\frac{\eta \left(T_f-t_1\right)}{{\sigma }_0}\right)\gg 1$

4) The stress increase rate η is proportional to the strain rate $\stackrel{\dot{}}{\epsilon }$ in [t₁, T_f].

$\eta =c\stackrel{\dot{}}{\epsilon }$ for $t>t_{\text{1}}$ until fracture .

To derive relation (5) under these conditions, we first formulate an approximate equation for the fracture time T_f. Using conditions 1 & 2, we express the contribution of σ(t) over [0, T_f] approximately as

$\frac{1}{t_0}{\displaystyle {\int }_0^{T_f}\exp \left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)\text{d}t}\approx \frac{1}{t_0}{\displaystyle {\int }_{t_1}^{T_f}\exp \left(\frac{w\left(t_1\right)+\eta \left(t-t_1\right)}{{\sigma }_0}\right)\text{d}t}$

$=\frac{1}{t_0}\exp \left(\frac{w\left(t_1\right)}{{\sigma }_0}\right)\frac{{\sigma }_0}{\eta }\left[\exp \left(\frac{\eta \left(T_f-t_1\right)}{{\sigma }_0}\right)-1\right]$

Condition 3 leads to approximation:

$\exp \left(\frac{\eta \left(T_f-t_1\right)}{{\sigma }_0}\right)-1\approx \exp \left(\frac{\eta \left(T_f-t_1\right)}{{\sigma }_0}\right)$

Thus, fracture time T_f satisfies approximately the equation

$\frac{{\sigma }_0}{t_0\eta }\exp \left(\frac{w\left(t_1\right)}{{\sigma }_0}\right)\exp \left(\frac{\eta \left(T_f-t_1\right)}{{\sigma }_0}\right)=1$

Solving for T_f, we obtain

$\left(T_f-t_1\right)=\frac{1}{\eta }\left({\sigma }_0\ln \left(\frac{t_0\eta }{{\sigma }_0}\right)-w\left(t_1\right)\right)$

The stress at fracture is

${\sigma }_f\left(\eta \right)={\left(w\left(t_1\right)+\eta \left(t-t_0\right)\right)|}_{t=T_f\left(\eta \right)}={\sigma }_0\ln \left(\frac{t_0\eta }{{\sigma }_0}\right)$

Using condition 4, $\eta =c\stackrel{\dot{}}{\epsilon }$, we write the stress at fracture as

${\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)={\sigma }_0\ln \left(\frac{t_{0}c\stackrel{\dot{}}{\epsilon }}{{\sigma }_0}\right)=\ln \left(10\right){\sigma }_0{\log }_{10}\left(\frac{t_{0}c}{{\sigma }_0}\right)+\ln \left(10\right){\sigma }_0{\log }_{10}(\epsilon \dot{})$

This is the same as equation (5), which was derived based on much stronger assumptions.

Appendix B: Fracture Tolerance vs Strain Rate in the Power Law Model

We consider the model in which the fracture time t_f(σ), as a function of applied static stress σ, follows a power law.

$t_f\left(\sigma \right)=t_0{\left(\frac{\sigma }{{\sigma }_0}\right)}^{-b}$

where

σ₀: a characteristic stress;

t₀: time to fracture at applied stress σ₀ (t₀ depending on the selection of σ₀).

With the power law model, the cumulative contribution toward fracture from the applied stress profile {σ(t)} over time interval [0, T] is

$\text{contributionof}\left\{\sigma \left(t\right),0<t<T\right\}=\frac{1}{t_0}{\displaystyle {\int }_0^T{\left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)}^b\text{d}t}$

We derive the relationship between the applied strain rate and the observed fracture stress when the stress and strain are locally proportional in a time interval leading to fracture. Specifically, we assume that the loading satisfies the conditions below:

1) The stress vs time profile is linear in a time interval leading to fracture. That is, σ(t) is linear for t > t₁ until fracture:

$\sigma \left(t\right)=\{\begin{array}{ll}w\left(t\right),\hfill & t<t_1\hfill \\ w\left(t_1\right)+\eta \left(t-t_1\right),\hfill & t>t_1\text{untilfracture}\hfill \end{array}$

2) Let T_f = the fracture time. The time interval [t₁, T_f] contains dominant contribution toward fracture. In other words, the contribution toward fracture from $\left\{\sigma \left(t\right),0<t<t_1\right\}$ is negligible:

$\frac{1}{t_0}{\displaystyle {\int }_0^{t_1}{\left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)}^b\text{d}t}\ll 1$

3) From t₁ to T_f, the contribution toward fracture per unit time increases significantly:

${\left(\frac{\sigma \left(T_f\right)}{\sigma \left(t_1\right)}\right)}^b={\left(1+\frac{\eta \left(T_f-t_1\right)}{w\left(t_1\right)}\right)}^b\gg 1$

4) In the time interval [t₁, T_f] leading to fracture, the stress increase rate η is locally proportional to the strain rate $\stackrel{\dot{}}{\epsilon }$ :

$\eta =c\stackrel{\dot{}}{\epsilon }$ for $t>t_{\text{1}}$ until fracture.

By the definition of fracture time T_f, the cumulative contribution of stress profile {σ(t)} over time interval [0, T_f] is 1. Using conditions 1 & 2, we write it as

$\begin{array}{c}1=\frac{1}{t_0}{\displaystyle {\int }_0^{T_f}{\left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)}^b\text{d}t}\approx \frac{1}{t_0}{\displaystyle {\int }_{t_1}^{T_f}{\left(\frac{w\left(t_1\right)+\eta \left(t-t_1\right)}{{\sigma }_0}\right)}^b\text{d}t}\\ =\frac{{\sigma }_0}{t_0\eta }{\left(\frac{w\left(t_1\right)}{{\sigma }_0}\right)}^{b+1}\frac{1}{b+1}\left[{\left(1+\frac{\eta \left(T_f-t_1\right)}{w\left(t_1\right)}\right)}^{b+1}-1\right]\end{array}$

Using condition 3, we further simplify it to

$1=\frac{1}{t_0}{\displaystyle {\int }_0^{T_f}{\left(\frac{\sigma \left(t\right)}{{\sigma }_0}\right)}^b\text{d}t}\approx \frac{1}{b+1}\frac{{\sigma }_0}{t_0\eta }{\left(\frac{w\left(t_1\right)+\eta \left(T_f-t_1\right)}{{\sigma }_0}\right)}^{b+1}$

Thus, the fracture time T_f satisfies approximately the equation

${\left(\frac{w\left(t_1\right)+\eta \left(T_f-t_1\right)}{{\sigma }_0}\right)}^{b+1}=\frac{\left(b+1\right)t_0\eta }{{\sigma }_0}$

Solving for T_f we obtain

$w\left(t_1\right)+\eta \left(T_f-t_1\right)={\sigma }_0{\left(\frac{\left(b+1\right)t_0\eta }{{\sigma }_0}\right)}^{\frac{1}{b+1}}$

The fracture tolerance, σ_f(η), is the stress at fracture.

${\sigma }_f\left(\eta \right)={\sigma \left(t\right)|}_{t=T_f}=w\left(t_1\right)+\eta \left(T_f-t_1\right)={\sigma }_0{\left(\frac{\left(b+1\right)t_0\eta }{{\sigma }_0}\right)}^{\frac{1}{b+1}}$

Now applying condition 4, $\eta =c\stackrel{\dot{}}{\epsilon }$, and taking the log on both sides, we obtain

${\log }_{10}{\sigma }_f\left(\stackrel{\dot{}}{\epsilon }\right)={\log }_{10}{\sigma }_0+\frac{1}{b+1}{\log }_{10}\left(\frac{c\left(b+1\right)t_0}{{\sigma }_0}\right)+\frac{1}{b+1}{\log }_{10}(\epsilon \dot{})$

This is Equation (9) in the main text.

Appendix C: Behavior of Q(Y)

For random variable Y, we consider quantity Q(Y) defined as

$Q\left(Y\right)\equiv 2{\left(\frac{E\left(Y\right)}{\text{std}\left(Y\right)}\right)}^2\left[\text{med}\left(\ln \left(Y\right)\right)-E\left(\ln \left(Y\right)\right)\right]$

When Y is a normal random variable with the mean significantly larger than the standard deviation to ensure Y being virtually positive, we write it as

$Y=Y_0+s_0\xi ,Y_0\gg s_0,\xi ~N\left(0,1\right)$

In this case, we derive Q ≈ 1.

$\text{med}\left(\ln \left(Y\right)\right)=\ln \left(\text{med}\left(Y\right)\right)=\ln (Y0)$

$\ln \left(Y\right)=\ln \left(Y_0+s_0\xi \right)\approx \ln \left(Y_0\right)+\left(s_0/Y_0\right)\xi -\frac{1}{2}{\left(s_0/Y_0\right)}^2{\xi }^2$

$E\left(\ln \left(Y\right)\right)\approx \ln \left(Y_0\right)-\frac{1}{2}{\left(s_0/Y_0\right)}^2$

$Q\left(Y\right)\equiv 2{\left(\frac{E\left(Y\right)}{\text{std}\left(Y\right)}\right)}^2\left[\text{med}\left(\ln \left(Y\right)\right)-E\left(\ln \left(Y\right)\right)\right]\approx 1$

When Y is a log-normal random variable, we have

$\text{ln}\left(Y\right)~\text{normaldistribution}$,

$\text{med}\left(\ln \left(Y\right)\right)=E\left(\ln (Y)\right)$

$Q\left(Y\right)\equiv 0$.

A log-normal distribution has a heavy tail. When the tail of Y is heavier than that of a log-normal, Q takes a negative value. For example, when Y is a log-lognormal, we have

$\text{ln}\left(Y\right)~\text{log-normaldistribution}~\text{exp}\left(\mu +\sigma N\left(0,\text{1}\right)\right)$

$\text{med}\left(\text{ln}\left(Y\right)\right)=\text{exp}(\mu )$

$E\left(\ln \left(Y\right)\right)=\exp \left(\mu +{\sigma }^2/2\right)$

$Q<0$

In summary, if Y is close to a normal distribution, we have Q ≈ 1; if Y is close to a log-normal distribution, we have Q ≈ 0; if the tail of Y is heavier than that of a log-normal distribution, we have Q < 0.

NOTES

¹Many tissue properties in the body are strain rate dependent, as are the dynamics of impact.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Shen, W., Niu, E., Webber, C., Huang, J. and Bykanova, L. (2012) Advanced Total Body Model (ATBM) Analyst’s Guide for Model Verification and Validation. L-3 Applied Technologies/Jaycor, Technical Report, No. J0939-10-389. https://mtec-sc.org/wp-content/uploads/2017/12/ATBM-Analysts-Guide-for-Model-VV-JNLW12-082v2.pdf
[2]	Wang, H., Burgei, W.A. and Zhou, H. (2018) Dose-Injury Relation as a Model for Uncertainty Propagation from Input Dose to Target Dose. American Journal of Operations Research, 8, 360-385. https://doi.org/10.4236/ajor.2018.85021
[3]	Pellettiere, J.A., Duma, S.M., Bass, C.R. and Crandall, J.R. (1998) Strength of the Female Upper Extremity. Proceedings NATO/RTO Specialist's Meeting: Models for Aircrew Safety Assessment: Uses, Limitations, Requirements, Dayton, OH, 26-28 October 1998, 7.1-7.12.
[4]	Duma, S.M., Schreiber, P.H., McMaster, J.D., Crandall, J.R., Bass, C.R. and Pilkey, W.D. (1999) Dynamic Injury Tolerances for Long Bones of the Female Upper Extremity. Journal of Anatomy, 194, 463-471. https://doi.org/10.1046/j.1469-7580.1999.19430463.x
[5]	Duma, S.M., Schreiber, P.H., McMaster, J.D., Crandall, J.R. and Bass, C.R. (2002) Fracture Tolerance of the Male Forearm: The Effect of Pronation versus Supination. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, 216, 649-654. https://doi.org/10.1177/095440700221600803
[6]	Sanborn, B., Gunnarsson, C.A., Foster, M., Moy, P. and Weerasooriya, T. (2014) Effect of Loading Rate and Orientation on the Compressive Response of Human Cortical Bone. Army Research Laboratory, Report No. ARL-TR-6907. https://www.arl.army.mil/arlreports/2014/ARL-TR-6907.pdf
[7]	Amitrano, D. and Helmstetter, A. (2006) Brittle Creep, Damage and Time to Failure in Rocks. Journal of Geophysical Research: Solid Earth, 111, B11201. https://doi.org/10.1029/2005JB004252
[8]	Carter, D.R. and Caler, W.E. (1985) A Cumulative Damage Model for Bone Fracture. Journal of Orthopaedic Research, 3, 84-90. https://doi.org/10.1002/jor.1100030110
[9]	Kress, T.A., Snider, J., Porta, D., Fuller, P.M., Wasserman, J.F. and Tucker, G.V. (1993) Human Femur Response to Impact Loading. Proceedings of the International IRCOBI Conference on the Biomechanics of Trauma, Eindhoven, 8-10 September 1993, 93-104.
[10]	Havaldar, R., Pilli, S.C. and Putti, B.B. (2014) Insights into the Effects of Tensile and Compressive Loadings on Human Femur Bone. Advanced Biomedical Research, 3, Article ID: 101. https://doi.org/10.4103/2277-9175.129375