+ μ 3 + σ 3 ε 3 (21)

where { ε 1 , ε 2 , ε 3 } are i.i.d. samples of N ( 0,1 ) . The binary injury outcome is governed by the sign of random variable

Z ( X ( t r u e ) , ω ) Z ( c ) = X ( e s t ) z ( c ) ( μ 1 + μ 2 + μ 3 ) + σ 1 2 + σ 2 2 + σ 3 2 ε (22)

At a given value of X ( e s t ) , random variable ( Z ( X ( t r u e ) , ω ) Z ( c ) ) has the same mathematical form as random variable ( Z ( x , ω ) z ( c ) ) in (11). As a result, the injury probability vs the estimated input dose has the expression

p | X ( e s t ) = x = Pr ( [ x z ( c ) ( μ 1 + μ 2 + μ 3 ) + σ 1 2 + σ 2 2 + σ 3 2 ε ] > 0 ) = 1 2 + 1 2 erf ( x z ( c ) ( μ 1 + μ 2 + μ 3 ) 2 σ 1 2 + σ 2 2 + σ 3 2 ) (23)

Injury function (23) has the same form as (12). Thus, p | X ( e s t ) = x is described by the normal distribution model with shape parameters ( D 50 , W ) given as follows.

p | X ( e s t ) = x = f normal ( x ; D 50 , W ) (24)

W = 2.584 σ 1 2 + σ 2 2 + σ 3 2

D 50 = z ( c ) + ( μ 1 + μ 2 + μ 3 )

In a well controlled lab setting, the true input dose X ( t r u e ) is measurable. For example, in experiments of male forearm fracture [12], a cylinder of specified mass is dropped from a specified height along a vertical track onto the PMHS forearm sample. Both the forearm sample and the cylinder are connected to accelerometers, allowing accurate measurements of the dynamic impactor load and the support loads. In addition, in situ strain gauges are used to record time series of strains at various locations during the loading. In this idealized setting, there is no measurement error in X ( t r u e ) . The injury probability as a function of the true input dose, p | X ( t r u e ) = x , can be determined from the observed binary injury outcomes vs measured values of X ( t r u e ) . Injury function p | X ( t r u e ) = x follows the normal distribution model with shape parameters ( D 50 ( 0 ) , W ( 0 ) ) given below.

p | X ( t r u e ) = x = f normal ( x ; D 50 ( 0 ) , W ( 0 ) ) (25)

W ( 0 ) = 2.584 σ 1 2 + σ 3 2

D 50 ( 0 ) = z ( c ) + ( μ 1 + μ 3 )

With this formulation, we can map back and forth between injury functions p | X ( t r u e ) = x and p | X ( e x t ) = x . We can also revise the injury function p | X ( e x t ) = x measured on one population to construct the injury function for a different population. We now discuss these two problems.

Problem 1:

Suppose we are given an injury model p | X ( t r u e ) = x , specified by shape parameters ( D 50 ( 0 ) , W ( 0 ) ) . The given injury function is for an idealized setting where the true input dose is directly measured. Our goal is to extend the given injury function p | X ( t r u e ) = x to predict the injury probability, p | X ( e x t ) = x , as a function of estimated input dose for the same population when the true input dose is not measurable.

Solution:

Injury function p | X ( t r u e ) = x is specified by shape parameters ( D 50 ( 0 ) , W ( 0 ) ) given in (25) while injury function p | X ( e x t ) = x is specified by shape parameters ( D 50 , W ) given in (24). Combining (25) with (24), we write ( D 50 , W ) as an update on ( D 50 ( 0 ) , W ( 0 ) ) .

W 2 = W 0 2 + ( 2.584 ) 2 σ 2 2 D 50 = D 50 ( 0 ) + μ 2 (26)

Problem 2:

Suppose we are given an injury model p | X ( e x t ) = x , specified by shape parameters ( D 50 , W ) . The given injury function is established based on measurements of a heterogeneous population, labeled population 1. Population 1 is characterized by uncertainties in the input dose estimation and in the critical threshold, as described in (20) and (21)

( X ( t r u e ) X ( e s t ) ) = N ( μ 2 , σ 2 2 )

Z ( c ) = z ( c ) + N ( μ 3 , σ 3 2 )

Now consider a different heterogeneous population, labeled population 2, with uncertainties described by

( X ( t r u e ) X ( e s t ) ) = N ( μ ˜ 2 , σ ˜ 2 2 )

Z ( c ) = z ( c ) + N ( μ ˜ 3 , σ ˜ 3 2 )

Here we assume that the propagation uncertainty from true input dose to target dose ( Z ( X ( t r u e ) , ω ) X ( t r u e ) ) = N ( μ 1 , σ 1 2 ) is statistically the same for the two populations. Our goal is to predict the injury function ( p | X ( e x t ) = x ) 2 for population 2 based on the given injury function p | X ( e x t ) = x for population 1.

Solution:

Injury function p | X ( e x t ) = x for population 1 is specified by shape parameters ( D 50 , W ) while injury function ( p | X ( e s t ) = x ) 2 for population 2 is specified by shape parameters ( D ˜ 50 , W ˜ ) . We write ( D ˜ 50 , W ˜ ) as an update on ( D 50 , W ) to take into account the differences in uncertainties between the two populations.

W ˜ 2 = W 2 + ( 2.584 ) 2 ( σ ˜ 2 2 σ 2 2 + σ ˜ 3 2 σ 3 2 ) D ˜ 50 = D 50 + ( μ ˜ 2 μ 2 ) + ( μ ˜ 3 μ 3 ) (27)

4. Dose-Injury Function for Target Doze of Log-Normal Distribution

For the discussion below, we adopt the normal-distribution model as the base formulation, switching away from the logistic model. There are several reasons behind the switching.

• The normal-distribution model is based on 1) viewing the binary injury outcome as completely determined by the target dose at the active site, 2) explaining the randomness in injury outcome as the consequence of uncertainty in dose propagation from input dose to target dose, and 3) modeling the dose propagation uncertainty as an additive Gaussian noise. This interpretation is both theoretically and operationally appealing.

• Mathematically, the injury function form of normal-distribution model is exactly invariant when additional normally distributed noise/uncertainty is incorporated into the model.

• We will study dose-injury models based on normally distributed intermediate variable. Mathematically, such an injury model is conveniently treated as a transformation of the normal-distribution model since the target doze is expressed as a function of the normally distributed intermediate variable.

• As we demonstrated in the previous section, the logistic model is practically equivalent to the normal-distribution model with the same shape parameters ( D 50 , W ) .

We first recall the function form of the normal-distribution model. In terms of internal variables ( σ , μ , z ( c ) ) , it is given by (12). In terms of shape parameters ( D 50 , W ) , it is expressed in (18). Geometric quantities D 50 , D 10 , D 90 and W of the injury function are related to internal variables ( σ , μ , z ( c ) ) as

D 50 = z ( c ) + μ D 10 = D 50 erf 1 ( 0.8 ) 2 σ D 90 = D 50 + erf 1 ( 0.8 ) 2 σ W D 90 D 10 = 2 2 erf 1 ( 0.8 ) σ (28)

Because of the symmetry of error function erf ( z ) , the normal distribution model (18) is symmetric around the median injury dose D 50 :

Symmetry : f normal ( D 50 + x ) 1 2 = 1 2 f normal ( D 50 x ) (29)

We now study a skewed injury function that breaks this symmetry. Consider the situation where the target dose Z ( x , ω ) has a log-normal distribution

Z ( x , ω ) = x e x p ( μ + σ ε )

Again ε ~ N ( 0,1 ) is a standard normal random variable. In this case, l n ( Z ( x , ω ) ) and l n ( x ) are simply related by an additive Gaussian noise.

l n ( Z ( x , ω ) ) = l n ( x ) μ + σ ε (30)

If we use l n ( x ) and l n ( Z ( x , ω ) ) to measure, respectively, the input dose and the target dose, then the injury probability vs l n ( x ) follows the same function form as (12) with ( x ; z ( c ) ) replaced by ( l n ( x ) ; l n ( z ( c ) ) ) :

p ( x ) = f normal ( l n ( x ) ; l n ( z ( c ) ) ) = 1 2 + 1 2 erf ( l n ( x ) l n ( z ( c ) ) μ 2 σ ) (31)

We examine the injury probability as a function of the original input dose x. The purpose is to investigate 1) under what condition the injury probability vs x can be approximated by the symmetric normal-distribution model, and 2) when the normal distribution approximation is invalid, what additional parameter we need to introduce to describe the injury function for the original input dose x.

Since the injury probability vs l n ( x ) follows the normal distribution model (12), we use results (28) for (12) to write out D 50 ( l n ) , D 10 ( l n ) and D 90 ( l n ) for quantity l n ( x ) .

D 50 ( ln ) = l n ( z ( c ) ) + μ D 10 ( ln ) = D 50 ( ln ) erf 1 ( 0.8 ) 2 σ D 90 ( ln ) = D 50 ( ln ) + erf 1 ( 0.8 ) 2 σ (32)

The corresponding D 50 , D 10 and D 90 for quantity x are

D 50 = e x p ( D 50 ( ln ) ) = z ( c ) e μ D 10 = e x p ( D 10 ( ln ) ) = D 50 e x p ( erf 1 ( 0.8 ) 2 σ ) D 90 = e x p ( D 90 ( ln ) ) = D 50 e x p ( erf 1 ( 0.8 ) 2 σ ) (33)

In this case, it is clear that ( D 90 D 50 ) > ( D 50 D 10 ) . The injury probability vs quantity x is not exactly symmetric around D 50 . We introduce a measure of skewness to represent the asymmetry of injury probability vs quantity x.

γ ln ( D 90 D 50 D 50 D 10 ) (34)

Specifically, γ defined above measures the skewness of interval [ D 10 , D 90 ] around D 50 .

• When γ = 0 , interval [ D 10 , D 90 ] is symmetric around D 50 .

• When γ > 0 , we have ( D 90 D 50 ) > ( D 50 D 10 ) , which implies that the upper half (above D 50 ) of injury function is flatter than the lower half (below D 50 ).

• When γ < 0 , we have ( D 90 D 50 ) < ( D 50 D 10 ) , and that the upper half of injury function is steeper than the lower half.

Skewness γ is an indicator of how well the injury function for x can be approximated by the symmetric normal distribution model. For a target dose of log-normal distribution, the skewness is γ = erf 1 ( 0.8 ) 2 σ . When σ is small, the skewness γ 0 , and the injury function is nearly symmetric around D 50 . When σ 2 > 0 , the skewness γ is positive, and in (31) the injury probability as a function of x is not symmetric. In this case, the injury function is characterized by three shape parameters: ( D 50 , W , γ ) .

D 50 = z ( c ) e μ W D 90 D 10 = D 50 2 s i n h ( erf 1 ( 0.8 ) 2 σ ) γ l n ( D 90 D 50 D 50 D 10 ) = erf 1 ( 0.8 ) 2 σ (35)

Notice that even though expressions of ( D 50 , W , γ ) in (35) contain three variables ( μ , σ , z ( c ) ) , two variables μ and z ( c ) appear only as a combination ( z ( c ) e μ ) in D 50 . Mathematically, the three shape parameters ( D 50 , W , γ ) are completely specified by ( D 50 , σ ) , and thus, have only two degrees of freedom. As a result, the three shape parameters ( D 50 , W , γ ) cannot be set independently of each other. For example, in (35) when γ is small, the width W will be small unless the median dose D 50 is large. Formulation (35), based on target dose of log-normal distribution (30), cannot accommodate any negative skewness ( γ < 0 ). It cannot even accommodate the simple symmetric case of γ = 0 with finite W > 0 and D 50 < + . We like to revise the formulation and construct an injury model in which the three shape parameters ( D 50 , W , γ ) can be set independently of each other.

5. A Dose-Injury Model with Skewness Based on a Normally Distributed Intermediate Variable

We construct a model that accommodates the median injury dose ( D 50 ), the width (W) and the skewness ( γ ) as 3 independent parameters. In previous section, we studied the formulation based on target dose of log-normal distribution, in which the skewness is always positive and the 3 shape parameters ( D 50 , W , γ ) are not independent of each other. A log-normal random variable can be viewed as the exponential of normal random variable. To accommodate negative skewness and to make ( D 50 , W , γ ) independent of each other, we extend the formulation to the case of target dose being a more general function of normal random variable.

We consider the situation where the dose propagation uncertainty is an additive Gaussian noise in quantity l n | x x 0 | with x 0 as a new tunable parameter. The target dose Z ( x , ω ) and the input dose x are related by

Z ( x , ω ) x 0 = ( x x 0 ) exp ( μ + σ ε )

In this setting, ( Z ( x , ω ) x 0 ) has the same sign as ( x x 0 ) . The domain of x is divided by x 0 into two regions: x > x 0 and x < x 0 . Only the region containing the critical threshold z ( c ) will be relevant for the injury model. The other region of x produces target dose Z ( x , ω ) always above or always below z ( c ) . For example, when x 0 > z ( c ) , only the region x < x 0 is relevant for the injury model; the region x > x 0 leads to target doze Z ( x , ω ) > x 0 > z ( c ) and thus, leads to an injury probability of 100%. We discuss separately the case of x 0 > z ( c ) and the case of x 0 < z ( c ) .

5.1. Case 1: x 0 < z (c)

In this case, the region x < x 0 yields target dose Z ( x , ω ) < x 0 < z ( c ) and an injury probability of 0%. We focus on the region x > x 0 , the relevant region for the injury model. The logarithm of shifted target dose l n ( Z ( x , ω ) x 0 ) and logarithm of shifted input dose l n ( x x 0 ) are related by an additive Gaussian noise.

ln ( Z ( x , ω ) x 0 ) = l n ( x x 0 ) μ + σ ε (36)

where ε ~ N ( 0,1 ) . We apply the shift x n e w = x o l d x 0 on all dose quantities (including D 50 and z ( c ) ). After the shift, problem (36) above is exactly the same as problem (30) in the previous section. It follows that the injury probability has the same function form as (12) with ( x ; z ( c ) ) replaced by ( l n ( x x 0 ) ; l n ( z ( c ) x 0 ) )

p ( x ) = f normal ( l n ( x x 0 ) ; l n ( z ( c ) x 0 ) ) for x 0 < z ( c ) = 1 2 + 1 2 erf ( l n ( x x 0 ) l n ( z ( c ) x 0 ) μ 2 σ ) (37)

Based on results (33) and (35), we write out ( D 50 , γ , W ) for injury function (37).

D 50 x 0 = ( z ( c ) x 0 ) e μ D 10 x 0 = ( D 50 x 0 ) e x p ( erf 1 ( 0.8 ) 2 σ ) D 90 x 0 = ( D 50 x 0 ) e x p ( + erf 1 ( 0.8 ) 2 σ ) W D 90 D 10 = ( D 50 x 0 ) 2 s i n h ( erf 1 ( 0.8 ) 2 σ ) γ l n ( D 90 D 50 D 50 D 10 ) = erf 1 ( 0.8 ) 2 σ (38)

Note that both z ( c ) and D 50 are on the right side of x 0 in the case of x 0 < z ( c ) . As we will see, D 50 and z ( c ) are always on the same side of x 0 . With Formulas (38) for the case of x 0 < z ( c ) , we can accommodate shape parameters ( D 50 , W , γ ) with positive skewness γ > 0 . Specifically, at any fixed μ , for each given set of ( D 50 , W , γ > 0 ) there is a unique corresponding set of ( x 0 , σ , z ( c ) ) .

σ = γ erf 1 ( 0.8 ) 2 ( D 50 x 0 ) = W 2 s i n h ( erf 1 ( 0.8 ) 2 σ ) , x 0 = D 50 ( D 50 x 0 ) . z ( c ) = x 0 + ( D 50 x 0 ) e μ (39)

This works for any positive skewness γ > 0 , corresponding to the situation where the injury probability has a flatter rise above the median injury dose D 50 than below it.

To accommodate negative skewness γ < 0 , however, we need x 0 > z ( c ) .

5.2. Case 2: x 0 > z (c)

In this case, we focus on the region x < x 0 since the region x > x 0 yields target dose Z ( x , ω ) > x 0 > z ( c ) and an injury probability of 100%. The target dose and input dose are related by

l n ( x 0 Z ( x , ω ) ) = l n ( x 0 x ) + μ + σ ε (40)

where ε ~ N ( 0,1 ) . Here we consider quantity l n ( x 0 Z ( x , ω ) ) with the negative sign because it is an increasing function of Z ( x , ω ) . Injury occurs when the target dose is above the critical threshold: Z ( x , ω ) > z ( c ) , which translates to

Injury probability = Pr ( l n ( x 0 Z ( x , ω ) ) > l n ( x 0 z ( c ) ) )

The injury probability has the expression

p ( x ) = f normal ( l n ( x 0 x ) ; l n ( x 0 z ( c ) ) ) for x 0 > z ( c ) = 1 2 + 1 2 erf ( l n ( x 0 x ) + l n ( x 0 z ( c ) ) + μ 2 σ ) (41)

Notice that (31) with quantities denoted by (' ) and (41) are connected by transformation

( x x 0 ) = 1 x , ( z ( c ) x 0 ) = 1 z ( c ) , μ = μ

We use results (33) and (35) for injury function (31) to write out ( D 50 , γ , W ) for (41).

D 50 x 0 = ( x 0 z ( c ) ) e μ < 0 D 10 x 0 = ( D 50 x 0 ) e x p ( erf 1 ( 0.8 ) 2 σ ) < 0 D 90 x 0 = ( D 50 x 0 ) e x p ( erf 1 ( 0.8 ) 2 σ ) < 0 W D 90 D 10 = ( D 50 x 0 ) 2 s i n h ( erf 1 ( 0.8 ) 2 σ ) > 0 γ l n ( D 90 D 50 D 50 D 10 ) = erf 1 ( 0.8 ) 2 σ < 0 (42)

In the case of x 0 > z ( c ) , both z ( c ) and D 50 are on the left side of x 0 . With Formulas (42) for the case of x 0 > z ( c ) , we can accommodate shape parameters ( D 50 , W , γ ) with negative skewness γ < 0 . Specifically, at any fixed μ , for each given set of ( D 50 , W , γ < 0 ) there is a unique corresponding set of ( x 0 , σ , z ( c ) ) .

σ = ( γ ) erf 1 ( 0.8 ) 2 ( D 50 x 0 ) = W 2 s i n h ( erf 1 ( 0.8 ) 2 σ ) , x 0 = D 50 ( D 50 x 0 ) . z ( c ) = x 0 + ( D 50 x 0 ) e μ (43)

This works for γ < 0 , which indicates that the injury probability has a steeper rise above the median injury dose D 50 than below it.

Next we combine the results of x 0 < z ( c ) and x 0 > z ( c ) to derive a unified formulation for accommodating shape parameters ( D 50 , W , γ ) regardless of the sign of γ .

5.3. A Unified Formulation for All Values of Skewness

In the previous sub-section, we studied models based on target dose of shifted log normal distribution with shift as a parameter. We now synthesize the results obtained to develop a unified formulation of injury function in which the 3 shape parameters ( D 50 , W , γ ) can be specified independently.

First, we show that at any fixed value of μ , there is one-to-one correspondence between ( x 0 , σ , z ( c ) ) and ( D 50 , W , γ ) . For any given set of shape parameters ( D 50 , W , γ ) regardless of the sign of γ , we combine results (39) and (43) to write out the corresponding ( x 0 , σ , z ( c ) ) .

σ = | γ | erf 1 ( 0.8 ) 2 x 0 = D 50 W 2 s i n h ( γ ) z ( c ) = x 0 + ( D 50 x 0 ) e μ (44)

Conversely, for any given set of ( x 0 , σ , z ( c ) ) , we combine results (38) and (42) to write out the corresponding shape parameters ( D 50 , W , γ ) .

D 50 = x 0 + ( z ( c ) x 0 ) e μ W = | D 50 x 0 | 2 s i n h ( erf 1 ( 0.8 ) 2 σ ) γ = sign ( z ( c ) x 0 ) erf 1 ( 0.8 ) 2 σ (45)

Again, D 50 and z ( c ) are always on the same side of x 0 . Next we combine (37) for x 0 < z ( c ) and (41) for x 0 > z ( c ) to write out a unified injury probability vs x.

p ( x ) = 1 2 + 1 2 erf ( sign ( z ( c ) x 0 ) 2 σ l n ( x x 0 ( z ( c ) x 0 ) e μ ) ) (46)

To specify the unified injury function in terms of shape parameters ( D 50 , W , γ ) , we express all quantities in (46) using only ( D 50 , W , γ ) and x.

sign ( z ( c ) x 0 ) 2 σ = erf 1 ( 0.8 ) γ

( z ( c ) x 0 ) e μ = ( D 50 x 0 ) = W 2 s i n h (γ)

x x 0 = ( x D 50 ) + ( D 50 x 0 ) = ( x D 50 ) + W 2 s i n h (γ)

With these expressions, we write the unified injury function as

p ( x ) = 1 2 + 1 2 erf ( erf 1 ( 0.8 ) γ l n ( 2 s i n h ( γ ) ( x D 50 ) W + 1 ) ) G ( x ; D 50 , W , γ ) (47)

In injury model (47), the 3 shape parameters ( D 50 , W , γ ) can be specified independently of each other. In particular, for small skewness γ 1 , expanding (47) in terms of γ reduces it to the symmetric normal-distribution model (18)

G ( x ; D 50 , W , γ ) = 1 2 + 1 2 erf ( 2 erf 1 ( 0.8 ) ( x D 50 ) W + O (γ) )

Figure 3 illustrates several injury functions of the form (47), respectively, for positive, zero and negative skewness. All injury functions shown have the same width W = 5 . In the left panel of Figure 3, injury functions are aligned at D 50 = 16 . This alignment demonstrates that for γ > 0 the left half of injury function is steeper than the right half; for γ < 0 the left half of injury function is flatter than the right half; and for γ = 0 the injury function is symmetric. In the right panel, injury functions are shifted to be aligned at D 10 = 13.5 and thus also aligned at D 90 = D 10 + W = 18.5 because they all have the same width

Figure 3. Injury functions with positive, zero and negative values of skewness. All injury functions have the same width W = 5 . Left panel: injury functions are aligned at the median injury dose D 50 = 16 . Right panel: injury functions are shifted to have the same interval [ D 10 , D 90 ] = [ 13.5,18.5 ] .

W = 5 . With D 10 = 13.5 fixed, the median dose D 50 varies with skewness γ from D 50 = 14.1 at γ = 2 , to D 50 = 16 at γ = 0 , and to D 50 = 17.9 at γ = 2 . The alignment of interval [ D 10 , D 90 ] highlights that as γ increases from negative to zero to positive, the injury function becomes more concave down.

6. Effect of Input Dose Uncertainty on the Injury Function with Skewness

We study the effect of input dose estimation uncertainty on the dose-injury function with skewness. We use the term “composite injury function” to denote the injury model after the input dose uncertainty has been incorporated into the model. In general, the composite injury function will be somewhat different from the 3-parameter function form (47) we derived in the previous section. We calculate the three shape parameters ( D 50 , W , γ ) of the composite injury function. Then we explore approximating the composite injury function using function form (47). We examine the difference between the composite injury function and model (47) with the same shape parameters ( D 50 , W , γ ) . If the approximation error is small, then the 3-parameter function form (47) is approximately invariant with respect to input dose uncertainty, and it serves as an adequate framework for accommodating uncertainty in estimating the input dose. Furthermore, framework (47) provides a mechanism of mapping the injury function for one particular dose propagation uncertainty to that for a different uncertainty. Using this mechanism, we can construct an injury model for a target population in application, based on measured injury data for a test population in experiments.

We start with a function of injury probability vs true input dose that is exactly of form (47) specified by 3 shape parameters ( D 50 , W , γ ) :

p 0 ( x ( t r u e ) ) = G ( x ( t r u e ) ; D 50 , W , γ )

We consider the situation where the true input dose x ( t r u e ) is not measurable. Instead, an estimated input dose, x, is obtained as an approximation for x ( t r u e ) . We assume

• the difference ( x ( t r u e ) x ) is a normal random variable, and

• the difference ( x ( t r u e ) x ) is independent of x.

We assess the injury probability as a function of the estimated input dose x. For each fixed value of x, the corresponding x ( t r u e ) is a normal random variable: x ( t r u e ) = x + μ + σ ε where ε ~ N ( 0,1 ) . The composite injury function, p σ ( x ) , representing the injury probability at estimated input dose x, is a Gaussian weighted average of p 0 ( x ( t r u e ) ) :

p σ ( x ) = E ( G ( x + μ + σ ε ; D 50 , W , γ ) ) = G ( x + μ + s ; D 50 , W , γ ) 1 2 π σ 2 exp ( s 2 2 σ 2 ) d s (48)

When injury function G ( ) has non-zero skewness, the Gaussian weighted average of G ( ) on the right hand side of (48) does not have a simple analytical expression. We use numerical integration to calculate the composite injury function p σ ( x ) and calculate its shape parameters ( D 50 ( σ ) , W ( σ ) , γ ( σ ) ) . We examine numerically if p σ ( x ) is still well described by function form (47) with p σ ( x ) ’s shape parameters ( D 50 ( σ ) , W ( σ ) , γ ( σ ) ) .

In our numerical study, p 0 ( x ( t r u e ) ) , the injury probability vs the true input dose before input dose uncertainty is incorporated, has function form (47) and is specified by shape parameters D 50 = 16 , W = 5 , and γ = 0.8 . We consider input dose uncertainty of normal distribution with μ = 0 (mean) and various values of σ (standard deviation). The composite injury function, p σ ( x ) , contains the effect of input dose uncertainty, showing injury probability vs estimated input dose x. Figure 4 examines the composite injury function p σ ( x ) for σ between 0 and 3.

The left panel of Figure 4 shows the injury probability vs the estimated input dose x, respectively, for σ = 0 , 1 , 2 , 3 . The most pronounced effect of input dose uncertainty is to spread out the injury function and increase the width. We examine the trend of shape parameters ( D 50 ( σ ) , W ( σ ) , γ ( σ ) ) when the input dose uncertainty σ is added and increased. The right panel shows ( D 50 ( σ ) , W ( σ ) , γ ( σ ) ) vs σ , of the composite injury function. As the input dose uncertainty σ increases, both the median injury dose D 50 ( σ ) and the width W ( σ ) increase monotonically, with W increasing more prominently than D 50 . At the same time, when σ increases, the asymmetry of injury function is smoothed out by the Gaussian noise and as a result, the skewness γ ( σ ) decreases. The change in median injury dose D 50 ( σ ) is attributed to the presence of skewness: the median injury dose increases (moves toward the right) when an injury function with positive skewness is smoothed out by a Gaussian noise. Conversely, the median injury dose decreases (moves toward the

Figure 4. Effect of the input dose uncertainty σ on the injury function with skewness. Left panel: composite injury functions for several values of σ . Right panel: shape parameters ( D 50 ( σ ) , W ( σ ) , γ ( σ ) ) vs σ of the composite injury function.

left) when an injury function with negative skewness is smoothed out. The movement of the median injury dose is caused by smoothing an asymmetric function (see Figure 3 for the general shape of injury functions with positive, zero, and negative skewness). For an injury function of zero skewness, D 50 ( σ ) is invariant with respect to σ when the injury function is smoothed out by a Gaussian noise.

Next we examine whether or not the composite injury functions for σ 2 > 0 shown in Figure 4 are still approximately described by model (47). Figure 5 compares the composite injury function and the approximation using function form (47) with shape parameters ( D 50 ( σ ) , W ( σ ) , γ ( σ ) ) of the composite injury function. The left panel of Figure 5 compares the composite injury function for σ = 1 and its approximation. The two functions are barely distinguishable from each other. To quantitatively examine the error of approximation, in the right panel we plot the difference between the composite injury function and its approximation. For all values of σ examined, the maximum error in approximation is less than 0.01 (1%). The results demonstrate that function form (47) specified by 3 independent shape parameters ( D 50 , W , γ ) is an adequate model for quantitatively describing general injury functions with skewness.

With the framework of function form (47) and mapping transformation (48), we can filter out the effect of input dose uncertainty in measured injury data. Suppose we are given a measured injury function, p σ 1 ( x ) , of form (47) for a particular population with input dose uncertainty σ 1 . We use transformation (48) to map it back to p 0 ( x ( t r u e ) ) , the injury function for the case of zero input dose uncertainty ( σ = 0 ). From there, we can apply the mapping transformation again to predict the injury model for another population with input dose uncertainty σ 2 . There is no simple analytical expression for the mapping

Figure 5. Approximation of the composite injury function p σ ( x ) using function form (47) with shape parameters ( D 50 ( σ ) , W ( σ ) , γ ( σ ) ) . Left panel: comparison of p σ ( x ) and its approximation for σ = 1 . Right panel: error of the approximation for several values of input dose uncertainty σ .

transformation. Both the forward and backward mappings need to be implemented numerically. The detailed numerical procedure will be discussed in a subsequent study.

7. Concluding Remarks

We considered injury models in the framework of dose propagation uncertainty. The mathematical formulation is based on that the binary injury outcome is completely determined by the target dose at the active site and the critical threshold. The randomness in the occurrence of injury at a given input dose is attributed to the dose propagation uncertainty from input dose to target dose. The normal distribution model describes the situation where the dose propagation uncertainty is normally distributed. We interpreted the widely used logistic model as a good approximation to the normal distribution model, and thus, interpreted it approximately as a consequence of normally distributed dose propagation uncertainty. In many applications, the input dose is not directly measurable. Instead, an estimated input dose is calculated via computer simulations from measured quantities using representative median parameter values of the general population. In many practical situations, injury models are constructed in the form of injury probability vs estimated input dose. The discrepancy between the estimated input dose and the true input dose can be viewed as an uncertainty in the input dose. With the interpretation of dose propagation uncertainty, the input dose uncertainty is conveniently incorporated into the injury model. The framework of dose propagation uncertainty provides a mechanism of extending an injury function established on a test population to predict the injury model for a different population in application. Both the logistic model and the normal distribution model are specified by two shape parameters: the median injury dose and the 10 - 90 percentile width. The mapping between the injury functions of two populations has a simple analytical form of updating the two shape parameters. Both the logistic model and the normal distribution model are symmetric around the median injury dose and have no skewness. To accommodate injury functions with skewness, we studied dose propagation uncertainties of shifted log normal distribution with shift as a parameter. Based on the shifted log normal model, we developed a function form for injury probability vs input dose that is specified by three shape parameters: median injury dose, the width, and the skewness. The proposed function form allows the three shape parameters to be set independent of each other. In particular, the proposed function form is capable of accommodating arbitrary skewness, positive or negative. In addition, we showed numerically that the proposed 3-parameter function form is approximately invariant with respect to additions or changes in input dose uncertainty. Therefore, the 3-parameter function form serves as a broad framework for modeling input dose uncertainty and modeling injury function skewness at the same time. This broad framework allows us to map injury function with skewness from a test population to a different population in applications.

Disclaimer and Acknowledgements

The authors thank C. Kramer and J. Swallow of Institute for Defense Analysis (IDA) for bringing the problem to their attention, and thank the Joint Non-Lethal Weapons Directorate of U.S. Department of Defense 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.

Cite this paper
Wang, H. , Burgei, W. 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. doi: 10.4236/ajor.2018.85021.
References
[1]   Vorst, M.V., Stuhmiller, J., Ho, K., Yoganandan, N. and Pintar, F. (2003) Statistically and Biomechanically Based Criterion for Impact-Induced Skull Fracture. Annual Proceedings/Association for the Advancement of Automotive Medicine, 47, 363-381.

[2]   Cazares, S.M., Finnin, M.S., Holzer, J.R., King, A.L. and Kramer, C.M. (2017) Significance of Rib Fractures Potentially Caused by Blunt Impact Non Lethal Weapons. Institute for Defense Analyses (IDA), Alexandria, VA.

[3]   Chan, P., Ho, K. and Ryan, A.F. (2016) Impulse Noise Injury Model. Military Medicine, 181, 59-69. https://doi.org/10.7205/MILMED-D-15-00139

[4]   Murphy, W.J., Khan, A. and Shaw, P.B. (2011) Analysis of Chinchilla Temporary and Permanent Threshold Shifts Following Impulsive Noise Exposure. Centers for Disease Control and Prevention, Cincinnati, OH.

[5]   Hampton, C.E. and Kleinberger, M. (2018) Material Models for the Human Torso Finite Element Model. US Army Research Laboratory (ARL), Aberdeen Proving Ground, MD.

[6]   Shen, W., Niu, E., Webber, C., Huang, J. and Bykanova, L. (2012) ATBM Analyst’s Guide for Model Verification and Validation. L-3 Applied Technologies, San Diego, CA.

[7]   Kramer, C. and Swallow, J. (2018) Error Propagation in RSI Estimates of Blunt Impact NLWs Using the Advanced Total Body Model. Institute for Defense Analyses (IDA), Alexandria, VA.

[8]   Cox, D.R. (1958) The Regression Analysis of Binary Sequences. Journal of the Royal Statistical Society. Series B, 20, 215-242.

[9]   Wang, H., Burgei, W.A. and Zhou, H. (2017) Interpreting Dose-Response Relation for Exposure to Multiple Sound Impulses in the Framework of Immunity. Health, 9, 1817-1842.
https://doi.org/10.4236/health.2017.913132

[10]   Sturdivan, L.M., Viano, D.C. and Champion, H.R. (2004) Analysis of Injury Criteria to Assess Chest and Abdominal Injury Risks in Blunt and Ballistic Impacts. The Journal of Trauma: Injury, Infection, and Critical Care, 56, 651-663.
https://doi.org/10.1097/01.TA.0000074108.36517.D4

[11]   Bliss, C.I. (1934) The Method of Probits. Science, 79, 38-39.
https://doi.org/10.1126/science.79.2037.38

[12]   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.

 
 
Top