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

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 elastic-ity-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


Introduction
Non-lethal blunt impact weapons have been widely used by law enforcement and military to incapacitate individuals while minimizing fatalities and collateral damage.  [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 1 and all of which is obtained on PMHS or animal sub- 1 Many tissue properties in the body are strain rate dependent, as are the dynamics of impact. 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 sult 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.

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

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 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.

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]: where D 50 = the median injury dose, at which injury probability p(D 50 ) = 50%; in general, D η = the dose value at which p(D η ) = η%; and W ≡ D 90 -D 10 = the 10 -90 percentile width of the injury function.

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: 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. • 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 (Zx) 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 (Zx) as an independent Gaussian random va- 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: The injury model is completely specified by two shape parameters: median injury dose D 50 and the width W ≡ D 90 -D 10 . It has the expression: where the shape parameters (D 50 , W) are related to (μ, s) as Both the logistic model and the normal distribution model are completely specified by shape parameters (D 50 , 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.

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 1 , 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.

Two Inference Formulations
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: We use the logistic function form ; , 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 50 , 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: 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 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:   Table 1 below compares the performance of the two methods in estimating the shape parameters (D 50 , W) from data. The data sets used in comparison study are generated by simulations using parameter values D 50 = 16 and W = 5. Note that mathematically we can shift and scale the dose-injury function to make D 50 = 0 and W = 1. The error study is not affected by the particular choice of (D 50 , 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:

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.

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 σ.
In dynamic loading experiments, the applied stress varies with time: σ(t). The cumulative contribution of applied stress profile 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: where η is the stress increase rate. We introduce time scale t * .  ( ) Solving for T f (η) we obtain In all meaningful dynamic loading experiments, we have That allows us to write the fracture time T f (η) approximately The fracture tolerance, σ f (η), is defined as the stress at fracture.
Let ε denote the strain increase rate (strain rate). Suppose the stress increase rate η is proportional to the strain rate ε until fracture, c η ε =  until fracture.
We write fracture tolerance σ f as a function of strain rate ε With the exponential model, the theoretically predicted fracture tolerance is a linear function of the logarithm of strain rate. In light of the linear relation of ( ) f σ ε  vs ( ) 10 log ε described in (5), we propose a loading rate adjusted dose for predicting fracture risk. 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.

Local Linear Relation between Stress Rate and Strain Rate
In the mathematical derivation above, we used the linear relation The theoretically predicted linear dependence of ( ) f σ ε  on ( ) 10 log ε , 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.

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 With the power law model, we show that the logarithm of fracture tolerance is 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 This loading rate adjusted dose is based on the power law model.

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 ε . Table 2 summarizes the test results reported in [6] for strain rates ranging from ε = 10 -3 /s to ε = 10 3 /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  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).
H. Y. Wang et al.  log log 1 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)σ 0 in (5), or coef- (9). For each of the compression directions, the values of ln(10)σ 0 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)σ 0 for the longitudinal compression is significantly higher than that for the trans-  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  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 3 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.

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.

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.
where E( ) denotes the mean and med( ) the median.
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,  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,  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.

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 50 , 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 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  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  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. Health   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.

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  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]. To connect real test data and idealized ATBM simulations in a proper doseinjury 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.
H. Y. Wang et al.

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 50 and width W. When applying the injury model to a general group of subjects, we will update the shape parameters (D 50 , 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 lognormal 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- 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.  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.

An Injury Model Based on Existing Data of Forearm Fracture
The forearm has a more complex structure than that of the upper arm. We use the fluctuation analysis approach developed in the previous section to study the forearm fracture test results in [4].

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 H. Y. Wang et al. 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.

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, 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.
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:  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.

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.

Separate Injury Models For Extension Fracture and Compression Fracture
In 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 50 ) is given by E(Z (c) ) and the width (W) is related to std(Z (c) ) by Extension fracture: x (E) ≡ ln(extension 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: is too low. In general, we expect the compression tolerance should be larger than the extension tolerance: . As a reasonable placeholder until a relevant estimate becomes available, we shall use 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.

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. We calculate the fracture risk corresponding to the strain distribution ( ) ; , where 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 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)    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 .

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 50 and the width W. The median injury dose D 50 is the mean of the non-observable ln(elasticity homogenized strain) at fracture. D 50 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 50 on age (year) and body mass (kg). 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

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.    shape parameters  of injury model  50  50 , , , , injury risk  50  for realization   50 , , , , 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 .
H. Y. Wang et al.

Options for Reducing the Number of ATBM Runs When Assessing the Injury Risk of a Population
Risk assessment calculation based on framework (13)    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:  Note that the approximation error in Monte Carlo integration is proportional to 1 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.

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) 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.

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 putational 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 σ ε σ ε σ 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 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 ( ) In this case, we derive Q ≈ 1.
med ln ln 1 std When Y is a log-normal random variable, we have 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 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.